A conservation law formulation of nonlinear elasticity in general relativity 
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We present a practical framework for ideal hyperelasticity in numerical relativity. For this purpose, 
we recast the formalism of Carter and Quintana as a set of Eulerian conservation laws in an arbitrary 
3+1 split of spacetime. The resulting equations are presented as an extension of the standard 
Valencia formalism for a perfect fluid, with additional terms in the stress-energy tensor, plus a set 
of kinematic conservation laws that evolve a configuration gradient ip A i. We prove that the equations 
can be made symmetric hyperbolic by suitable constraint additions, at least in a neighbourhood of 
the unsheared state. We discuss the Newtonian limit of our formalism and its relation to a second 
formalism also used in Newtonian elasticity. We validate our framework by numerically solving a 
set of Riemann problems in Minkowski spacetime, as well as Newtonian ones from the literature. 
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Neutron stars are believed to form a crystalline outer 
crust as they age and cool, but retain a fluid (probably 
superfluid) core [1] . A mathematical framework for weak 
solutions of general relativistic elasticity is likely to be 
indispensable for the modelling of neutron star crusts in 
at least two scenarios: star quakes and binary mergers. 

Pulsars are observed to spin down at a regular rate, 
losing angular momentum through gravitational and/or 
electromagnetic radiation. Occasionally the rotation 
spins up suddenly. One model suggests that such a 
"glitch" occurs when the elastic crust breaks and the in- 
ertial moment of the star decreases suddenly as a conse- 
quence (e.g. [2]). It has also been suggested [3 that star- 
quakes are the cause of soft gamma repeaters (SGRs). 
Quasi-periodic oscillations in the tails of giant flares in 
SGRs have been suggested [4] to provide direct obser- 
vational evidence for crust oscillation modes, although 



the modelling of neutron star oscillations even in per- 
turbation theory is complicated by the coupling between 
the crust, the fluid core and a strong magnetic field (see 
e.g. [5]). A correct model would of course have to be non- 
linear. Finally we note that strong shocks also arise when 
two old neutron stars in a binary system merge. The de- 
tailed dynamics and features, such as the breaking strain 
(see [6 ), of the crust, will determine when and where 
the crust melts and breaks. This will in turn impact on 
the post-merger dynamics, such as the time taken by the 
remnant to collapse to a black hole (see e.g. [7]). 

For all these scenarios, models must therefore comprise 
an elastic crust, a fluid core, and a magnetic field perme- 
ating both. As a step towards such models, we present 
here a formulation of (hyper) elastic matter in general rel- 
ativity in the form of conservation laws amenable to solu- 
tion by high-resolution shock-capturing (HRSC) numer- 
ical methods. These conservation laws are the union of 
the usual stress-energy conservation (dynamics), and a 
set of conservation laws for a deformation tensor (kine- 
matics). 

The kinematic equations are essentially the same in 
Newtonian and relativist ic physics, but the literature on 
weak solutions of Newtonian elasticity uses Cartesian 
tensor notation, which obscures the geometric nature of 
the theory. In Sec. [TTJ we derive these equations care- 
fully, using the language of differential geometry. Follow- 
ing Carter and Quintana [8 , we begin with a map from 
spacetime to a 3-dimensional matter space. The main 
object we calculate with is its derivative ty A i. As a par- 
tial derivative, this is subject to integrability conditions. 
Under a 3+1 split these become evolution equations and 
constraints, of a purely kinematic nature, both of which 
can be written as conservation laws. We show that their 
physical significance is to allow discontinuities in the den- 
sity and kinks but forbid discontinuities in the crystal 
axes and particle world lines. 

The other, dynamical, half of the problem consists in 
finding the stress-energy tensor from ip A i and an equa- 
tion of state. We do this in Sec. |III[ following Karlovini 
and Samuelsson [9]. In particular, demanding covariance 
on both spacetime and matter space restricts the pos- 
sible dynamics. For isotropic matter, the equation of 
state can relate only two deformation scalars, besides the 
number density, internal energy, and entropy. Similarly, 
the stress-energy tensor depends on the equation of state 
through two generalised forces. 



In Sec. IV, we prove that our kinematic and dynamical 
equations together form a first-order system of evolution 
equations that, by constraint addition, can be made hy- 
perbolic if the constraints (12) are obeyed or not. This 



property is crucial for the stability of numerical solu- 
tions in which the constraints are left to evolve freely, 
and hence finite difference error generically generates 
constraint violations. We use the methods of Beig and 
Schmidt [10 , who proved symmetric hyperbolicity of an 
inequivalent first-order system (one which would not be 
appropriate for modelling weak solutions). 



In order to make contact with existing work on ideal 
fluid dynamics and magnetohydrodynamics in general 
relativity, in Sec. [V] we present our dynamical equations 
as a generalisation of the Valencia [11 formulation of hy- 
drodynamics. We give an algorithm for the conversion 
between conserved and primitive variables. 

As a first test of our formalism, we present numerical 



time evolutions of Riemann problems in Sec VI The vari- 
ables are three-dimensional, and the grid is either one- 
dimensional, or two-dimensional with the Riemann prob- 
lem at an angle to the grid. We compare the relativistic 
code in the Newtonian limit with an explicitly Newtonian 
code, and both with published Newtonian Riemann prob- 
lems [I2j [13] . We also compare against exact Riemann 
solutions in the relativistic regime in Minkowski space- 
time. We compare the Eulerian and mixed formalisms, 
and evolutions where the number density is either read 
off from the deformation tensor, or evolved separately. 



In Sec. VII we summarize the results of the paper, and 
discuss the work remaining to apply these methods to 
full 3+1 nonlinear simulations. 

We collect relevant formulas from the standard 3+1 
split of spacetime in Appendix [A) and relevant defini- 
tions of hyperbolicity in Appendix|B| One of two existing 
Newtonian formalisms [HI [15] is essentially the Newto- 
nian limit of our formalism. We derive the Newtonian 
limit in Appendix [C] In Appendix [D] we derive the equa- 
tions of an alternative Newtonian formalism [16-18 , and 
prove that the two have the same weak solutions. 

The remaining Appendixes contain auxiliary material 
on our numerical method and our numerical tests: Ap- 
pendix [E] proposes a general framework for discrete con- 
straint preservation (similar to "constrained transport" 
for MHD), and Appendix |f| presents our implementation 
of Riemann tests on a 2-dimensional grid. Appendix |H| 
describes the equations of state we use. Appendix [I] sum- 
marizes how we construct exact solutions for specific Rie- 
mann problems and Appendix [j] the initial data for our 
Riemann tests used here. 

We have attempted as far as possible compatibility 
with the notation of [9], [10] and [TT] , Throughout this 
paper, tensor indices are assumed to be in a generic local 
coordinate basis. Partial derivatives in this basis are indi- 
cated by commas. Indices a, 6, c, • • • = 0, 1, 2, 3 are space- 
time indices, z,j, &, ••• = 1,2,3 are spatial indices on 
x° = t = const, hyper surf aces, and A,B,C, — • = 1, 2, 3 
are matter space indices on a 3-dimensional matter space 
X3. In Appendix [Dj the indices a, /3, 7, • • • = 0, 1, 2, 3 are 
matter space indices on an extended matter space X 4 . 
In Sees. |III| and \LV\ a, j3 = 1,2 are used to label elas- 
tic forces. In Appendix [Bj a, j3 label the variables of a 
generic hyperbolic system. For all these indices a sum- 
mation convention applies. 

In order to take determinants of 2-index objects which 
are not (1,1) tensors, we introduce the non-tensorial to- 
tally antisymmetric symbol £, which is defined to be 
^0123 = 1, etc. With the exception of the objects £, 
throughout this paper, all objects transform as tensors 



of the type indicated by their free indices, unless we indi- 
cate otherwise by a suffix: for example, the determinant 
of the spacetime metric in coordinates x a will be denoted 
by -g x - 



II. KINEMATICS 

A. The configuration gradient and its 3+1 split 

In the relativistic framework of [HJ [9], the matter con- 
figuration is encoded in a map from 4-dimensional space- 
time to 3-dimensional matter space 



X- M 4 ^X 3 , 



(1) 



or in local coordinates x a on spacetime and £, A on matter 
space, 



1°^^ = X A (x a )- 



(2) 



For simplicity of notation we denote the derivative d\ of 
X by a new symbol -0, 



ip: M 4 ^TX 3 (g}T*M 4 , 



x a i->^ 



dx a ' 



(3) 



(4) 



For time evolutions, we introduce a time- foliation of 
the spacetime, so that we have 



X : R x M 3 -»• X 3 , 


is) 


(*,**) -^ 


(6) 


with derivatives 




* < ■= w' ^ : = ir- 


(7) 



Following [19 , we shall call \ A the configuration and 
both ifj A a and ifj A i the configuration gradient. 

The matter space coordinates ^ A label particles and 
must therefore be constant along particle world lines, so 
that 



U a ^ A a = 0, 



(8) 



where the 4- velocity u a is tangential to the matter world 
lines. Parameterising the 4-velocity in the standard way 
as 

u a = (u\u i ) = a- 1 W(l,v i ) (9) 

(see Appendix |A| for more details), we have 

$\ = ~vH A i. (10) 

The configuration gradient ifj A i is also used as the pri- 
mary variable in the Newtonian framework of [14] [15] [21] 
(denoted there by g). This framework is the Newtonian 



limit of our relativistic one. In Appendix [C] we derive 
the Newtonian limit of our framework. Other Newtonian 
papers [17] [18] use the the 3x3 matrix inverse of i\) A i, 
which we shall denote by F l a, as the primary variable 
(denoted there by_F). We review this alternative frame- 
work in Appendix [d] In the Newtonian literature, F % a is 
commonly called the (Lagrangian) deformation gradient, 
and ip A i the inverse deformation gradient. From a geo- 
metric point of view, however, these objects on their own 
carry no information about what one might intuitively 
call a deformation. 



B. Kinematic equations and hyperbolicity fix 

From the definition of ifj A a as a partial derivative, we 
have the integrability conditions 



C ab'-=i> [a,b] 

In a 3+1 split, these become 



0. 



E' 



2C A it = V\ t + (v^ A j) 



0. 



(11) 



(12) 
(13) 



The constraints ( 12 ) are conserved by the evolution equa- 



tions (13). Note that these equations are already in con- 



servation law form: more explicitly, 



E 



A 



v\ f + (*>'vV?) k = °- 



(14) 
(15) 



Instead of E A i = as an evolution equation for i\) A i, 
we shall in fact use 



?A 



E A { := 2aW- l u a C A ia = 2^ M +2v>tl> A lid] = 0. (16) 
This can be written as a balance law obtained from the 



conservation law (13) by adding a source term that is 



proportional to the constraint ([12]), namely 

^ A l ,t+(v j ^ A J ) l = 2v^ A M . 



(17) 



Note that this cannot be written in pure conservation law 
form. 

In handwaving anticipation of the hyperbolicity analy- 
sis presented in Sec. IV, we point out in passing that (17) 
can be written as an advection equation for ip A i with a 
source term that is of lower order in ip A i, namely 



ip A i + + v j ip A i, 



-lp A jV\ 



(18) 



For given v z , this is strongly hyperbolic in i\) A i, whereas 
(13) is only weakly hyperbolic. 



C. Kinematic jump conditions 



The geometric meaning of the integrability conditions 



(11) is that the particle world lines and the instanta- 
neous crystal lines (i.e. lines of constant £ A ) all mesh 
up into a four-dimensional grid. In particular, the world 
lines and crystal lines are continuous. The weak form 
of these equations must therefore keep them continu- 
ous while allowing them to kink, thus forbidding dislo- 
cations and fractures. To stress their purely kinematic 
nature, we shall discuss them without invoking a metric 
on spacetime or matter space. We ignore the source term 
in the evolution equations (17) in deriving the Rankine- 



Hugoniot conditions, because it has no effect on physical 
solutions, which obey the constraints. 

Consider a surface of discontinuity in space (from now 
on called a shock for briefness). Let ni be a covector nor- 
mal to the shock (uniquely defined up to an overall fac- 
tor). Let s l be the shock velocity vector (defined, in the 
absence of a metric, only up to the addition of a vector 
tangential to the shock), and let s := s l rii (which inher- 
its the arbitrary factor in rii but not the arbitary vector 
in s l ) be the normal shock speed. The jump (Rankine- 
Hugoniot) conditions arising from ( 14 ) and (p~5J) are then 



v>V^] 



-s [v> 



ni =0, 
1 1 + \v j ^ A M n k = 0. 



3i . 



(19) 
(20) 



We want to decompose these conditions into parts nor- 
mal and parallel to the shock. Let n % be a vector that 
obeys n l rii — 1. ri 1 is therefore uniquely defined up to 
the factor in n^, and the addition of an arbitrary vector 
tangent to the shock. Define the tensor 



?* 



n n 



j- 



(21) 



It is the projection operator into the tangent plane of the 
shock in the sense that 



■■ rii = 0, 



0, 



ik * =ir fc - 



(22) 



Split into normal and tangential components defined by 
v n := v l m and v" z :=|| 2 jV^ : the jump conditions can now 
be compactly written as 



w 



o, 



[lP A „(v n -s)]+Tp A lli [v\\ i }=0. 



(23) 

(24) 



The first of these guarantees the continuity of crystal lines 
(£ A lines) across the shock, or the absence of "surgery 
across the shock" , as illustrated in Figs, [I] and [2] The sec- 
ond guarantees the conservation of particles as they cross 
the shock. Consider the special case where v^ % is contin- 
uous. Then, in the rest frame of the shock, [ip A n v n ] = 0. 
This is a pure "density" shock of the type familiar from 
fluid dynamics, and is illustrated in Fig. [3J Conversely, 
consider the case where ^ A n is continuous. Then, again 
in the rest frame of the shock, ip A n [v n ] + ^yJ'O"*] = 0. 



FIG. 1. A discontinuity of the type illustrated here is not 
allowed by the jump conditions (it would require "surgery" on 
the material). For simplicity and without loss of generality, 
we choose space and matter space coordinates in this and 
the next three figures so that the shock is along the y axis 
and ip A i = Sf in the left state. This type of surgery is then 



forbidden by [ip 



= 0. 




FIG. 2. This type of surgery is forbidden by [^ 
coordinate choices are as in Fig. [I] 



0. The 



This is a pure travelling kink, set up by a discontinuity 
in the tangential velocity, as illustrated in Fig. |4j 

Fluids allow for a contact discontinuity where the tan- 
gential velocity jumps. This is replaced by the travelling 
kink in elastic matter. (The only contact discontinuity 
that survives is the one where the entropy jumps.) This 
holds even in the limit where the dynamics goes to the 
fluid limit (the stiffness goes to zero and the stress-energy 
tensor becomes that of a fluid), and so the fluid limit is 
singular. 



D. Matter space metric and particle number 
current 



The minimal geometric structure on matter space is a 
volume form uabc whose integration over a volume in 
matter space gives the number of particles in that part 
of matter space. In addition, at least a conformal metric 
is required to define angles on matter space, which can 
then be compared with angles on spacetime to define 
deformations. But together these two structures define a 



Formally, tensor fields on matter space could be defined 
as tensors on spacetime whose Lie derivative along u a 
and contractions with u a all vanish, and this is indeed 
the approach of |5], and partly of [9]. However, equiv- 
alently the components kAB(x C ( xd )) can be considered 
as scalars on spacetime that are constant along particle 
world lines, so that 



u a k AB n = 0, 



(28) 



or in coordinates 



FIG. 3. A pure density shock, with v y continuous, shown 
in the rest frame of the shock. This is essentially a one- 
dimensional phenomenon. "Density" and velocity on the left 
and right (shown as arrows) are related (in these coordinates) 
by [^ X x v x ] = 0. 




FIG. 4. A pure travelling kink, shown in the rest frame of 
the left state. For simplicity, we have assumed ip x x = 1 to be 
continuous, which implies that v x = is continuous and the 
"volume density" is continous. However, the "line density" 
along the Y crystal axis is discontinuous. The shock speed 
and shear speed (shown as arrows) in these coordinates are 
related by s(ip Y x ) R = (v y ) R . 



full Riemannian metric fc^e- ("Distances" are measured 
in particles, not meters). Therefore we now assume that 
kAB is defined and uabc is compatible with it. In matter 
coordinates ^ A this means that 



riABC = yk^dABC, 



where 



5 ABC 5 DEF k AD k BE k CF 



(25) 



(26) 



is the usual determinant. The suffix £ is a reminder that 
it is not a scalar on matter space but depends on the £ A 
coordinates. 

We use kAB as an example to discuss the "evolution" of 
tensors on matter space. Matter space itself has no time, 
but as we are using a Eulerian framework, we effectively 
consider kAB(x C ' (x d )) as a function on spacetime. The 
push-forward of kAB to a tensor k a b on spacetime obeys 



£ukab = 0, u a k ab = u b k ab = 0. 



(27) 



kAB,t + v l k A B, 



0. 



(29) 



Numerically, we prefer to work with fc^e, which has fewer 
components and a simpler evolution equation than k ao . 
Following [8, 9 , we consider the push- forward of tlabc 
to a 3-form n a 6 C on spacetime 



riabc := i) A a ^ B h i> C \n A BC- 



(30) 



Spacetime also has a volume form e a 5 CC ^, compatible with 
a Lorentzian metric g^. In arbitrary coordinates, 



^abcd 



Jx Odbcd") 



■.= -^5 abcd S e ^ h g ae g bf g cg g dh . 

(31) 
(We have defined g x as positive for ease of notation) . We 
then define the particle number current 



j a := ^e abcd n bcd . 



This is timelike, and conserved, 
V a j° 



abcdyj 



o, 



(32) 



(33) 



where V a is the covariant derivative compatible with 
g^. The right-hand side vanishes because it is the push- 
forward of ri[BCD,A]i which must vanish as it is a 4-form 
on a 3-dimensional space. We split j a into a matter 4- 
velocity and a particle density 



j =: nu 



where u a is normalised as 



u u„ 



(34) 



(35) 



(and hence n = —j a j a )- In coordinates, using (A2) and 
( |A12| ), V a j a = becomes 



{s/TxWn) t + (y^Wntf) = 0. 



(36) 



Conversely, we can relate the partic le den sity and cur- 
rent via n = —u a j a , and substituting (A13) into this and 



using (A3), we obtain 
1 



3! 



W- X e ijk 



Tlijk 






(37) 



Here n^ are the space components of the 4-dimensional 
3- form riab c in the adapted coordinates (£, x % \ and ip x £ is 
the determinant 

^« := ^S^SABC^i^j^k- (38) 

We now show explicitly that V ' a j a = is a linear com- 
bination of the evolution equations (18) for i\) A i, that is, 



the kinematic evolution equations with the hyperbolicity 
fix. Contracting (18) with F l a-, the matrix inverse of 
ip A i, and using the matrix identity 5(lm/j x ^) = F l A^ A i-> 
we obtain 



(ln^Xt + fi^ln^Xi + fi* 



0. 



(39) 



Working from the other end, we insert (37) and ([9| into 
( 34 ) and use ( A2 ) to obtain 



k -0(^). 



/9x 
Hence V ' a j a = is equivalent to 

But with the advection equation 

(**),* + v i (kz) ti =0, 
which follows from (|29|), this is equivalent to (39) 



(40) 

(41) 
(42) 



III. RELATIVISTIC DYNAMICS 
A. Action and stress-energy tensor 

We begin with the matter action 

S:= fe(g ab ^ A a ,k AB ,...,s)g 1 x / 2 d i x, (43) 



where the dots stand for any other tensors on matter 
space and s is the entropy per rest mass (a scalar on mat- 
ter space). Varying for now only the metric, the standard 
definition of the stress-energy tensor T a 5, 



1 



SS=:- J T ab 5g ab g"J l ^x 



evaluates to 



lab ~ 2 dg^b~ e9ab - 



(44) 



(45) 



We define a projector into the tangent space normal to 
the 4- velocity, 



h a b '= u a u h + g ah . 



(46) 



hab should not be confused with the proje ctor ^ a b into 
the t = const hypersurfaces defined in (A7). 



We can now write 



where 



Tab = eU a U b +Pa6, 



9 de h 

Pab '' = dg^~ a6 ' 



(47) 



(48) 



which is by definition symmetric. 

We define the pull-back of the spacetime metric to mat- 
ter space, 



g AB :=r a ^ B b 9 ab - 



(49) 



We define qab as its matrix inverse. We therefore now 
have two Riemannian metrics on matter space, namely 
gAB and Jzab • As a matter of convention and terminol- 
ogy, we will refer to Jvab (only) as the matter space met- 
ric, but we will later implicitly move matter space indices 
(only) with gAB and g AB . Note that in this convention 
k AB := g AC g BD kBD, and that this is not the matrix in- 
verse of kAB- (We note in passing that the Newtonian 
limit ip A i ip B •7*- 7 ' of g AB is commonly called the Finger 
tensor in the Newtonian literature. The Newtonian lit- 
erature implicitly assumes that kAB and 7^ are flat and 
given in Cartesian coordinates and moves indices implic- 
itly. Moreover, some expressions can only be made sense 
of if F % a and ifj A i are also used implicitly to convert be- 
tween space and matter space indices.) 

As a further illustration of these conventions, the quan- 
tity 



^A a := ^ B b9 



ab 



9AB 



(50) 



is the inverse of ip A a (which is not a square matrix, and 
so has no matrix inverse) in the sense that 



^ A ai>B a 
^ A a^A b 



S A 



K b . 



(51) 
(52) 



(The first of these follows directly from the definition of 
gAB as the matrix inverse of g AB . The second can be 
shown by verifying that the right-hand side is normal to 
u a and u 01 and obeys h a b h c = h a c .) 

From covariance in both spacetime and matter space, 
we must have 



e(^ a ,g ab ) = e(g^), 



(53) 



as this is the only way the spacetime indices on ip A a and 
g ab can be contracted. (A more formal proof is given in 
[10].) Hence 



de 



de dg 



AB 



de 



Q Q ab Q g AB Q g a 



dg 



AB 



r a r 



b- 



(54) 



Hence p a bU a = 0, and so u a hb c T ab = 0. This means that 
there is no energy flux relative to the matter. In this 
sense we are dealing with ideal (non-dissipative) elastic 
matter. p a b is called the pressure tensor (for a perfect 



fluid, p ab = ph a01 where p is the pressure), and we now 
see that the Lagrangian e in the action ( [43] ) evaluates (for 
solutions to the Euler-Lagrange equations) to the total 
energy density (in the rest frame of the matter). 
We next note that 



n 2 = ^n abc n abc 



= ^g ad g be g cf n abc n def . (55) 



From its relation to the matter space volume form (30), 
n aoc is independent of g ab in the sense that it is con- 
structed only from tiabc and ty A a . Hence, taking a 
derivative of (55), 



dn 



dff 



1 



ab 



nh ab 



(56) 



where in the partial derivative n is considered as a func- 
tion of g ab , ifj A a and the matter tensors, as well as s. 
Then, defining e by 



=:rc(l + e), 



we have 



Pab = 2n 



de 



dg 



ab' 



(57) 



(58) 



with the same definition of the partial derivative. ([9] 
and [10] define e = ne. Here we take the rest mass out of 
the energy density to agree with the usual definition of e 
in relativistic hydrodynamics as the internal energy per 
rest mass.) Similarly to (54), we can write (58) as 



Pab = nr A B^ A a^ B b, 



where we have defined 



tab '•— 2 



de 



dg 



AB ' 



(59) 



(60) 



(The Newtonian limit of tab is commonly called the sec- 
ond Piaola-Kirchhoff tensor in the Newtonian literature, 
modulo the implicit assumptions mentioned above.) 



B. Isotropic matter 

We now specialise to the case that the specific inter- 
nal energy e depends on g ab , i/j A a , s and a single mat- 
ter tensor, the metric &ab- (Modelling matter with an 
anisotropic crystal structure would require e to depend 
on additional tensor fields on matter space, such as a pre- 
ferred frame.) e and hence e should transform as a scalar 
both on spacetime and on matter space. We therefore 
need to find all double scalars that can be made from 
g AB and Uab- 

From (49), we see that g AB transforms as a (2,0)-tensor 
on matter space and as a scalar on spacetime. With this 
in mind we define 



g AC k BC = g ac iP A a ip C ck 



Be- 



rn 



This transforms as a scalar on spacetime and as a (1,1) 
tensor on matter space. Hence its eigenvalues transform 
as scalars on matter space. They are the required double 
scalars. (We note that [9 work with the (l,l)-tensor on 
spacetime k a b = g aCr ^ B b^ C c^bc instead. This has the 
same eigenvalues as k A b plus one zero eigenvalue.) 

We split the matrix k A B into its determinant k and a 
unit determinant matrix n A B , 



n A - 

V B •" 



k-^k" 



(62) 



and note that the determinant is related to the particle 
density by 

1 

3f 



k:~- 



J A Bc5 DEF k A D k B E k c F 



5 B A 5 E B 5 F c] k A D k B E k c F 

g [AlD k AD g^ E k BE g lc]F k CF 
1 

3!< 



i AD g BE g CF n A Bcn DE F 



= n 2 , (63) 

where the first equality is the usual definition of the de- 
terminant of a matrix, the second reminds us that for 
a (l,l)-tensor this is actually a scalar, the third is the 
definition of k A b-, the fourth follows from the fact that 



tiabc is the volume form of kAB-, and the last one is (55) 
pulled back to matter space. 

We can now consider the specific internal energy e as 
a function of n, t] A b and s. In fact, it can depend on 
r] A B only through its scalar invariants, of which there 
are precisely two independent ones. Hence 

e(k A B , s) = e(k, V A B , s) = t{n, I 1 , I 2 , s), (64) 

where n = fc 1 / 2 as just shown and we have denned 



k-'/ s g AB k AB , 



(65) 



I 1 

I 2 ■= V A bV B a = k-^g AB 9 CD k AC k BD . (66) 
With gAB defined as the matrix inverse of g AB we have 
dk 



AB 



and hence 



dg 



dn 

AB 



kgAB, 



1 



dg 



-ng A B- 



We find 



TAB = -gAB + 2(/i7T^ B + /2^Ab)^ 



where 



p 

A,2 



^AB 



kab 



n 



2^_ 

dn 
de 



a/1,2' 

di 1 

dg AB ! 

dP_ 

~AB ' 



VAB ~ -gABl 1 , 



dg 



ZiVACV B - ^gABl )• 



(67) 

(68) 

(69) 

(70) 
(71) 
(72) 
(73) 



Substituting (69) into (59), we see that 



IV. HYPERBOLICITY 



Pab =ph a b +7Ta6, 



(74) 



with the first term the stress tensor of a perfect fluid and 
the second term representing the anisotropic stress, 

7Ta6 =^ A a^ B b^AB, (75) 

where 

7T AB := 2n(f 1 7T 1 AB + H^ab)- (76) 

Hence 7r a 5 is a tracefree spatial tensor in the sense that 

7T ab U a = 0, h ab 7T ab = 0. (77) 

Moreover, 7r a 5 vanishes if e depends only on n and s, 
which is the fluid limit. 

We also note that with the temperature defined by 



de 
da' 



(78) 



the first law of thermodynamics on a per particle basis 
can be written as 

de = Tds-pd(-j + /i dl 1 + f 2 dl 2 , (79) 

so /i } 2 are "generalised forces" in the thermodynamical 
sense. 



C. The unsheared state 

Elastic matter at a given density n has an unsheared 
state that minimises e at fixed n, but one cannot assume 
that there exists a relaxed state that minimises e abso- 
lutely, including under variation of n. This is because at 
sufficiently low pressure, and hence n, the matter may be 
in a fluid rather than solid state [9 . 

It is intuitively clear that the unsheared state corre- 
sponds to r] A B = 5 A b- In fact, we see from ( 72|73 ) 
that 7r a fr vanishes for all values of ip A a if and only if 



7j b = &B- This means that t]ab is the matrix inverse 
of g 



AB 



or 



Hence 



Vab = 9ab- 



kAB = n 2/s g AB 



(80) 



(81) 



in the unsheared state. It is natural to assume that mat- 
ter freezes in the unsheared state. Hence we set kAB to 



(81) at the moment of freezing, and advect it via (29) 
afterwards. Note that gAB is the pull-back of /i a 6, which 
even in special relativity is not flat, so in general kAB 
will not be flat, except in the Newtonian limit where 
hab — lab is flat and even then only if n takes a constant 
value at freezing. 



A. Overview 

For smooth solutions, it is natural to consider the rel- 
ativistic elasticity equations as a system of second-order 
PDEs in the variables x A ( x<1 )' I n order to show exis- 
tence and uniqueness of solutions, Beig and Schmidt [10] 
have introduced an explicit reduction to first order of 
these equations, and have shown that the reduction is a 
first-order symmetric hyperbolic system, at least in the 
unsheared state. 

The reduction of any second-order system to first-order 
hyperbolic form is complicated by the fact that the reduc- 
tion creates definition constraints on the auxiliary vari- 
ables (here, ip A [ij] = 0), which can be added to the 
evolution equations to change their principal part and 
hence their hyperbolicity properties. (We note in this 
context that in [20] a definition of symmetric hyperbol- 
icity for a second-order system has been given as the 
existence of a symmetric hyperbolic reduction to first, 
together with a necessary and sufficient criterion for this 
reduction to exist, which is purely algebraic in terms of 
the principal symbol of the second-order system. Hence 
if well-posedness of the second-order system is the only 
concern, constructing an explicit first-order reduction is 
unnecessary.) 

We have a different reason for constructing an explicit 
first-order reduction: we want to construct a numerical 
scheme that can accurately reproduce weak solutions of 
the relativistic elasticity equations. As for weak solutions 
of fluid mechanics, the standard way of doing this is to 
construct HRSC numerical schemes for the equations in 
an appropriate first-order balance law form. 

In this section we will show that the kinematic evolu- 



tion equations (17), together with the dynamical evolu- 
tion equations VbT ab = (constraints), form a symmet- 
ric hyperbolic system of evolution equations for ip A a , 
or equivalently ip A i and v % , if the constraints (12) are 
obeyed or not. We also show that (17), together with 



just V&T a6 = 0, as used in the Newtonian formalisms 



[1214141 121] , is strongly hyperbolic but not symmetric hy- 
perbolic. 

For completeness, relevant standard definitions of hy- 
perbolicity are summarised in Appendix El 



B. The second-order system 

Roughly speaking, the first-order equations for i(j A a 
must be the second-order equations for \ A \ replacing 
X A ,ab by ip A a> 5 and adding multiples of the constraint 
^ A [a,6] to the right-hand sides. We therefore derive the 
second-order equations first, following [10]. In particular, 
this will allow us to establish the standard connection be- 
tween the matter evolution equations and stress-energy 
conservation. 



Hence, in this subsection we consider g ab and \ A as 
the independent variables. We consider ifj A a = x A ,a as 
a derived object, and we consider fc^, any other matter 
space tensors, and s, as fixed tensor fields on matter space 
that are not varied in the following. The action is 



■■=Je(g ab , 



X A ,X A ,a,k A B,---,s)^d 4 x. (82) 



After integration by parts, and neglecting the boundary 
terms, its variation is 



SS 



1 



Tab Sg aD + S A W 






(83) 



where the stress-energy tensor is given as before by (45), 
and the Euler-Lagrange equations are 



Sa- 



de 
dXA 



1 

[9x 



/Qx 



de 



d(x A ,a) 



(84) 



Note that these are second-order differential equations 
for \ A • 

Variations generated by an infinitesimal change of co- 
ordinates x a — >• x a + ( a on the spacetime take the form 

SX A = £<X A = C¥ A C , (85) 

5g ab = C c g ab = 2V {a ( b K (86) 

The action must be invariant under such changes, and 
hence after another integration by parts 



V b T ab = ip A a £ A . 



(87) 



Hence stress-energy conservation holds if and only if the 
elastic matter field equations hold. 

£ a = has only three independent components, while 
V 6 T a 5 = has four. However, i^V^X^ = is equivalent 
to 



n 



(e-^n)+^ a6 V a ^ = 0, 



where a dot denotes u a V ' a . But this is just the first law 
(79), evaluated along a particle worldline, for smooth so- 



lutions, so that 5 = 0. Hence it is an identity if the stress- 
energy tensor is thermodynamically consistent with the 
equation of state. 

The matter equations in their second-order form can 
be written as 



£A = M a0 ABX U M-GA = Q 



where 



M ab 



d 2 e 



AB '■- 



d^ A a d^ B b 



(89) 



(90) 



are the coefficients of the principal part and Ga com- 
prises all lower-order terms. From its definition, 



M ab AB = M ba BA . 



(91) 



We shall see that M ab ab as defined by ( |90| ) is not sym- 
metric in ab alone, even though M^ab does not con- 
tribute to ([89]). 



C. The principal symbol 

We shall write the principal symbol more explicitly in 
terms of the shear and the equation of state. From (53), 



M ab AB - 4- 



d 2 e 



QgACQgBD 

where we have defined 



r 



^Ca^Db + 2 



4> A b g ab - 



de 



dg 



^,ab 



AB- 



With 



-« a « 6 + fTV B , 



(92) 



(93) 



(94) 



this can be split into parts parallel and normal to the 
4-velocity as 



M ab AB = -VA B u a u b + Uacbd^V, (95) 



where 



V-ab '■= 2 
Uacbd '■= 4 



de 
d 2 e 



(96) 



de 



dg AC dg BD ■ 2 ^ab9cd- (97) 



Note that there are no cross terms, that is Uah^M^ ab = 
0. 

We now evaluate the symbols \iab and Uacbd further. 
With (57), using (68), we can rewrite 



Vab = nTAB + egAB, (98) 

Uacbd = n(gAC^BD + Qbd^ac + t A b9cd 

+t A cbd) + 1eg A [ C gD]B, (99) 



where tab was defined above in (60), and analogously 
we have defined 



tabcd '•= 4 



d 2 e 



ABpi n CD ' 



dg AB dg' 



(100) 



Using the chain rule, we now express tab and tabcd 
as a sum of terms, each of which is a product of a matter 
scalar (such as e, p, c 2 s etc.) and a tensor that depends 
only on the deformation. We can rewite the expression 



(69) for tab more compactly as 



TAB = -QAB + l/a^ABi 

n 



(101) 



where a = 1,2 labels the shear scalars, and we use a 
summation convention over a. With the same notation, 
we can write 

P , / 2 P\ 
tabcd = —*—9a(c9d)b + \c s J 9ab9cd 

+2n (gABfna^CD + QCdJuol^'ab) 
^fa^AB^CD + ^ABCD, (102) 
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where 



02, 



2 = dp , ere 



fa/3 ■- 



and 



(103) 



"ABCD ■- dg AB dg CD 



( l 1 A rl 

I t:9a(c9d)b + q9ab9cd I J 



-^(i1ab9cd + Vcd9ab), (104) 



2 7-2 



2 _ d*f 

"ABCD ■- dg AB dg CD 



/ 2 4 \ r2 

= ( z9a(c9d)b + -^9ab9cd I -« 

;{r)AEn E ' b9cd + riCEV E d9ab) 



+2Va(cVd)b- 



(105) 



D. The unsheared state 

The principal symbol simplifies considerably in the un- 
sheared state, denoted by a circle, where 



I a = 3, 



(106) 
(107) 



^ABCD = 9A(C9D)B - T.9AB9CD, (108) 



11 ABC D — ^ABCDi 

and therefore 



(109) 



where 



riT AB =pgAB, (HO) 

riTABCD = 2rg A (c9D)B + Q9ab9cd, (HI) 



r:=-2> + 2n(/i+4/ 2 ), (112) 

g:=nc2-2>-^n(/i+4/ 2 ). (113) 

Note that in the unsheared state only the combination 
/i + 4/ 2 appears. 
We finally obtain 

Pab = (p J re)g A B, (114) 

Uacbd = (2p + g + e)g A c9BD + (p + r)g A BgcD 

+ (r - e)g A D9BC- (115) 

(This reduces to Eq. (4.16) of [10] in the special case 
p = e = 0.) 



E. First-order systems 

Any first-order reduction of the second-order system 
must have the form 



with 



where 



M ab AB^ B b„ 



ab 



■G, 







^06 



M a0 AB := M a0 AB + D a0 AB 



V AB — JJ l ] AB 



(116) 
(117) 
(118) 



governs constraint addition. 

In particular, VbT ab = should give us the dynamical 
part of the equations, but as we shall see, in order to 
achieve symmetric hyperbolicity of the entire system, we 
will have to add constraints to these equations. It turns 
out that adding constraints only to the "spatial" part of 
M ab AB is sufficient. Hence, we consider the evolution 
equations 



E a . = y^ab + A ACBD ^a^JJ^B [bc] = Q? (ng) 

where 

^ACBD = —J^ADBC (120) 

parameterises a family of constraint additions. To write 
VbT ab in terms of ip A a , we start from 



Tab = 2 ^ Aa ^ Bb 



eg 



ab 



(121) 



Keeping only the principal part in the matter variables, 
that is terms of the form ip B b,c-, we find after some cal- 
culation that 



V h T ab = f tp Aa M cb AB 



4^^y ]a 



dg 



AB 



^ B 6,c + 1.0. 

(122) 



Substituting (|94|) into (|122|), we find that 

E a = U Aa l 
where 



X M' 



cb 



AB' 



-2u 



> A ^ A V^\ c + l.o., (123) 



M ab AB 



M ab 



AB + (^ACDB ~ 2g AlC fJL D ] B )^ °T ■ 

(124) 
The modification of M ab ab can be pulled back to a mod- 
ification of Uacbd, namely 

Uacbd •= Uacbd + ^acdb — ^gA[c^D\B (125) 

Splitting E a into its parts parallel and normal to the 
4- velocity, we have 

u a E a = -2fjLAB4> Ab u c rl> B lbtC] , 



^ A aE Q 



M c 



AB^ b,c 



G A . 



(126) 

(127) 



But if \±ab is invertible, u a E a = is equivalent to the 
kinematic evolution equations with hyperbolicity fix (17). 
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F. Symmetric hyperbolicity 

The definition of symmetric hyperbolicity for a gen- 
eral system of first-order evolution equations is reviewed 
in Appendix [B] Roughly speaking, the principal symbol 
must be symmetric, and its time component must be pos- 
itive definite. We begin with the first condition. 

Symmetry The principal symbol of neither (116) nor 
of (123) has the correct index structure P a f3 c in the com- 
posite index defined by w a := i^ A a - We follow the ap- 
proach of Beig and Schmidt [10 . Define 

2 M b]a BA (128) 



W ab AB c := u a M cb AB - 2^' 



and consider the system of first-order equations 
S a A := W ab AB c ^j B b , c - G A u a = 0, 
A A := -u c g A BX B ,c = 0, 

where \ A and ^ A a are now considered as 
dent varia bles, and u A is defined by the i(j a 



(129) 

(130) 

indepen- 
through 
(30J32J34). It is clear that each solution \ A of the second- 
order system generates a solution ip A a = x A ,a of this 
first-order system. Beig and Schmidt [10] prove the con- 
verse, that a solution i(j A i of the first-order system obey- 
ing i/) A [ij] = gives rise to a second-order solution \ A • 

The principal symbol of this system has now the cor- 
rect index structure. It is easy to see that it is symmetric, 
in the sense that 



W a ° AB C = W° a BA 



ba 



-u c g A B = ~u c g B A, (131) 



if and only if M ab ab has the symmetry (91). To achieve 
this, we set 

^ACBD = J^ADBC + ^9A[C^D]Bi (132) 

where Aadbc has the symmetries 

^ACBD = ^BDAC = ~ ^ADBC, (133) 

and will be determined when we consider positivity of 

the principal symbol below. 

We now verify that the system ( 129|130 ) is eq uivalent 
to our evolution equations. From (95) and (124), we find 

u a M ab 



AB =U jlAB (134) 

and hence 

u a W ab AB C = -M cb AB , (135) 

^CaW ab AB c = -2^ d K c ^Uacbd. (136) 

If Uacbd is invertible as a matrix with composite indices 
AC and BD, we finally have the decomposition 

u a 8 a A = o M cb AB^ B b , c = G A , 

l^ E a £ a A =0 & "V[a,6]=0, 

< 



A A 



u a X C , 



o, 



We see that (137) and (138) are the same as rtl26 



in our formalism. Finally, ( 139) is equivalent top8j) (plus 




similar equations for any other matter tensors), as 

U a k A B,a = k A B,CU a X C ,a (140) 

by the chain rule. 



Positive definiteness The second condition for sym- 
metric hyperbolicity is the existence of a timelike covec- 
tor t a which makes the quadratic form (energy norm) 



E := t c {W ab AB C m A a m B h - u c g AB l A l 



AjB\ 



(141) 



positive definite, called a subcharacteristic vector. (The 
formal arguments m A a and l A of this quadratic form are 
in the tangent bundle of the phase space, and can be 
thought of as perturbations of ip A a and x A about a back- 
ground solution.) Decomposing m A a uniquely as 



m 



=: a A u a + a AC i/j C a, 



(142) 



and choosing t a = u a , we have 

E = fi AB a A a B + U A CB D a AC a BD + g AB l A l B . (143) 

Hence u a is a subcharacteristic vector if \±ab and Uacbd 
are positive definite. (Note that the n the y are also invert- 
ible, as we assumed earlier.) From (114), fi AB is positive 
definite in the unsheared state if p + e > 0. 

It remains to look at the positive definiteness of 
Uacbd- We choose 

^acbd = 2(d - e - p)g A [c9D]B, (144) 

or equivalent ly 

^acbd = ^nf a g A [c^D]B + 2d 9A[c9D]B (145) 

for the total constraint addition, where d is a parameter 
to be determined now. For simplicity, we look again at 
the unstrained case. Uniquely decomposing a AB as 



a AB = ^AB + ^AB 



;9 



AB 



(146) 



where the first term is antisymmetric and the second 
symmetric and tracefree, we find 



f T AC„BD 

Uacbd® ol 



nc 2 s + — j k 2 + duj AB uj A B 
[An(h+Af s )-d]K AB K AB - 



(147) 

We now see that this quadratic form is positive definite, 
and hence our evolution equations are symmetric hyper- 
bolic, for < d < 4n(/i + 4/ 2 ) (assuming that c 2 s > 0). 
Hence, adding some constraints to V&T a6 = is neces- 
sary for symmetric hyperbolicity, for example with the 
mid-range value of d = 2n(fi + 4/2). 

Tracing all the definitions back, we can write this par- 
ticular constraint addition as 

E a = V b T ab +h ac ^ Db [2nf a 7T% B + 4n(/i + 4f 2 )g DB ] $* c] 

(148) 
Looking back, the first term in the square brackets makes 
the principal part of the second-order system symmetric, 
and the second makes it positive definite. 
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G. Characteristics of the first-order system 



As reviewed in Appendix [Bj k a is a characteristic cov- 

ector of the first-order equations with characteristic vari- 

m B b 



able w^ = m B h if 



W ab AB c m B b k c 



0. 



(149) 



Once again we decompose m 5 in the form (142). Fix- 



ing an irrelevant overall factor, we parameterise the wave 
number k a as 



k a — XlL a e a 



(150) 



where e a = ip A a £A is a unit covector on spacetime nor- 
mal to u a and ca the corresponding unit (with respect 
to g AB ) covector on matter space. As reviewed in Ap- 
pendix |Bj A is then the physical velocity of the mode rel - 
ative to the matter. Using the decomposition ( 135|136 ), 
(149) is equivalent to the pair 



B e D 



Xa 



BD\ 



Uacbd {ol 
Xiiabol 3 + U A CBDe c a BD 



(151) 

(152) 



Moreover, symmetri c hy perbolicity implies that Uacbd 
is invert ible and so (151) is equivalent to 



a B e D 



Xa 



BD 



0. 



(153) 



Eq. (153) has two classes of solutions. One class obeys 



A = 0, 



a 



0, 



with a BD restricted by (152) to obey 



Uacbdz ol 



c^bd 



0. 



(154) 



(155) 



These modes travel at zero speed relative to the mat- 
ter. As (155) represents 3 equations for 9 components of 



a BD , there are 6 such modes. They can be parameterised 



(156) 



explicitly as 






a BD = {V -l )A CBD VAWc 


w c e c = 


0. 


The other class obeys 






A^0, a BD = - 


-X~ 1 a B e D , 




or equivalently 






m B h = A _1 a B fc 6 . 





(157) 



(158) 



Hence these modes are physical, obeying the constraints. 
Furthermore, all physical modes are of this form, which 
indicates that the modes in the class (154) are all un- 



physical. Substituting (158) into ( |152[ ), we find 

(-X 2 liab + UACB D e c e D )a B = 0, (159) 

or equivalently 

(-X 2 fi AB + U A CB D e c e D )a B = (160) 



(constraint addition drops out). This can be written as 

A AB a B = 0, A AB := M ab AB k a k h . (161) 

But, as reviewed in Appendix [Bj this is precisely the con- 
dition for k a to be a characteristic covector of the second- 
order system with characteristic variable a A . Hence the 
physical modes of the first-order system correspond one- 
to-one to the modes of the second-order system. There 
are 6 of these, forming 3 pairs with speeds d=A relative to 
the matter. 

H. Characteristics of the second-order system 



We now look at the solutions of ( 161 ) in more detail. 



The general expression for A^^ is quite long, and so we 
begin our analysis with the unsheared state. We find 



A A b = (-AX 2 + B)g AB + Ce A e B , 



(162) 



where 



A = e+p, B=p + r, C = 2p + r + q. (163) 

We can now read off the characteristic covectors and 
characteristic variables by inspection. Transversal waves 
have eigenvectors obeying ol b e B = (and so have two 
polarisations travelling with the same velocity), and 
Aab& B = then gives 



A : 



2 _ B _ 2/i + 8/2 ,2 

rT • A' 



A 1 



T- 



(164) 



Longitudinal waves have eigenvectors a B oc e B and 
Aab(* b = gives 



X 2_B±C c\ 4 2 _ 2 



(165) 



Taking the Newtonian limit of these characteristic 
speeds, we can identify the shear modulus \i and the bulk 
modulus K as 



K 



n(2/i 
nc 2 a . 



■8/2), 



(166) 
(167) 



(These expressions hold in units where the speed of light 
c is one, and where n is the rest mass density, rather than 
the particle number density. Otherwise they have to be 
multiplied by c 2 and the particle mass.) 

In the general, sheared, case the matter space tensor 
Aab is constructed from qab^ 9 AB \ Va B and e^. It 
would therefore be natural to decompose ca (and a B ) 
into eigenvectors of t]a B , which are automatically also 
eigenvectors of qa b = $a B - This can be done trivially 
by assuming that the index A labels that basis, so that 
t]a B is diagonal. The result is of the form A(A) = A 2 A 2 + 
Ao (dropping the indices on Aab)- We have solved the 
resulting cubic equation for A 2 by computer algebra but 
the result is too complicated to be illuminating. 
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Furthermore, for numerical purposes we are interested 
in the coordinate speeds A = dx z /dt in the x\ i = 1, 2, 3 
direction of some Eulerian coordinate system, that is, in 
characteristic covectors of the form 



k a = X(di) a - {dx l ) a . 



(168) 



This gives a characteristic equation of the form 

Aa = (A 2 A 2 + AAi + A )a = 0. (169) 

The resulting A are not related to the characteristic 
speeds relative to the matter 4- velocity in a simple way, 
because of the appearance of the Loren tz fac tor W in the 
relativistic velocity addition. To solve (169) numerically, 
we use a standard linear algebra package to find the right 
eigenvectors (Aa, a) T and eigenvalues A of the matrix 



-(A 2 )" 1 A 1 -(A 2 )~ 1 A 1 
I 



(170) 



I. Strong hyperbolicity 

In the Newtonian literature, the evolution equations 
for ifj A i are taken to be (17) (with constraint addition), 
but no constraints are added to VbT ab = 0. Our re- 
sults above show that the first-order system is then defi- 
nitely not symmetric hyperbolic, as the term Uacbd in 
the principal symbol is then not positive definite even 
in the unsheared state, and is not symmetric in general 
(although it is symmetric in the unsheared state). How- 
ever, the first-order system is still strongly hyperbolic if 
it admits a complete set of characteristic variables. We 
have just done the calculation in Sec. IV G[ and need to 
see only what changes if we cannot assume that Uacbd 
is symm etric and positive definite. 



(151) is no longer equivalent to (153) because Uacbd 



may not have an inverse, but solutions of (153) are solu- 



tions of (151). Furthermore, (155) only admits a larger 
solution space if Uacbd does not have maximal rank, so 
it will still have 6 solutions, even if they can no longer be 
explicitly parameterised by ( |156| ). Hence the 6 unphys- 
ical modes still exist. The 3 physical modes are com- 
pletely unaffected by con strain t addition (as one would 
expect) because we solve (160) to find them. 

Hence we have proved that the first-order system con- 
sisting of (17) and V a T ab = is strongly hyperbolic (but 



not symmetric hyperbolic) as long as the second-order 
system is strongly hyperbolic. 



V. STRESS-ENERGY CONSERVATION IN 3+1 
FORM 

A. Conservation laws 

The energy and momentum conservation laws in gen- 
eral relativity are the spacelike and timelike components 



of stress-energy conservation 

V a T ab = 0. (171) 

In [11], this is decomposed into four balance laws as 

V a [T a \d j ) b \=T ab V (a {d J ) b) , (172) 

V a {-T ab n b )=-T ab V (a n b) , (173) 

where n a is the unit normal on the t = const surfaces. 
The right-hand sides can be seen as a failure of stress- 
energy conservation to split into separate proper conser- 
vation laws for the energy and each momentum compo- 
nent, due to the failure of the thee spatial coordinate 
basis vectors (dj) a and the timelike unit normal vector 
n a to be Killing. (The choice of the fo ur ba sis ve ctors 
is merely conventional). In a 3+1 split, (172) and (173) 
become 



(a 2 ^T x T 00 ) t + (a 2 ^f x T M ) 



(174) 
(175) 



We now restrict to the elastic matter stress-energy ten- 



sor 



nab 



(e+p)u a u"+pg ab + n 



ab 



(176) 



To insert this into ( |174|) a nd (175K we need certain com- 
ponents of 7r ab . From (|75|) and (10) we have 



7T00 = 7TijV Z V J , 



TTOi = —KijV' 



n0 



(177) 



Using the 3+1 decomp osition o f the metric, the compo- 
nents that we need in ( 174|175 ) are 



n 00 = a- 2 -K ij v i v j , 


(178) 


ir°i = a~ 1 ir ij v i , 


(179) 


ir i j = h ik -a- 1 p i v k )ir kj , 


(180) 


TT 0i = (a-y-? - a - 2 p i v i )v k ir jk . 


(181) 


= 0, we have 




V l V 3 'Kij = ■f^TTij =: 7T. 


(182) 



From g ab 7r ab 



In numerical hydrodynamics, the conservation laws for 
the stress-energy are closed by the equation of state to- 
gether with the explicit particle number conservation law 



V a (mx a )=0. 



(183) 



As we have shown in Sec. |II Pi this evolution equation 
for n is equivalent to that for xjj^i (with the hyperbolicity 
fix) even when the constraints are not obeyed. Therefore, 
where we write n below, either value could be used with- 
out changing the hyperbolicity. However, we shall test 
this numerically by obtaining a value n^ from ip A i and a 
value no from D, and using either the one or the other. 
The conservation laws ( 172|173|183 ) can be written in 
the form 



(V7^0 t ~^~ { a \/lx^ 1 ) ■ = solirc e- terms. 



(184) 
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(Note the explicit insertion of y^ and a, which is only 
a convention - we follow [11 ). The conserved variables 
U = (D, Si,r) are related to the primitive variables via 



D = anu° = nW, 


(185) 


Si = aT°i = nhW 2 Vi + qtt * 




= nhW 2 Vi + nijV^ 


(186) 


r = a 2 T 00 -D = nhW 2 -p + a V - 


-D 


= n(hW 2 - W)-(p-ir), 


(187) 



where we have defined the standard specific enthalpy 

h := 1 +e+ -■ (188) 



p 

n 



The corresponding fluxes are 

F(D) 1 = nu l =na~ 1 Wv\ 
FiSjY = T l j = nhW 2 ^ 1 ^ +p5 i j + 
F(tY = aT 0i - F(D) 

= n(hW 2 - WOorV -hpoT 1 /^ 



j> 



n(W 2 - W)^- 1 ^ 

+(p - 7r)a -1 ^ + J ij 7r jk v 



(189) 
(190) 



(191) 



Following [11 , we have subtracted the rest energy from 
the total energy density in order to obtain the usual New- 
tonian energy conservation law in the Newtonian limit. 



B. Conversion of conserved to primitive variables 

In any numerical scheme, we frequently need to calcu- 
late the primitive variables n, v % and e from the related 
conserved variables D, Si and r, and hence further vari- 
ables such as p and 7r^ that appear in the fluxes. We 
assume that we have an equation of state that relates e, 
n, s and 7 a , and that we can use this to find the gener- 
alised forces p and f a . We also assume that g^ is evolved 
using the Einstein equations, that Jzab is advected using 
(29) and that ifj A i is evolved using the balance law (17). 
(Note that for our purposes ipf is both a conserved and 
a primitive variable.) 



The obvious difficulty is that to calculate TTij from ttab, 
we need v % while 7r^- is needed to extract v % from Sj. 
We therefore need to proceed iteratively. We guess the 
primitive variables 



From 



p 



r + D + (p - tt) = nhW 2 =: Z, 

(Si-TTikV^Sj-TTjiV 1 )^ 1 = Z 2 V 2 



(192) 



(193) 
(194) 



we obtain Z and v and hence W from ( A15 ), followed by 



riB from (185) and n^ from (37) (we compute both, but 
choose one value to use as n in the remainder of the cal- 



culation), v l from (186) and h from (187). We now have 



will not be consistent with the equation of state. We 
therefore now recompute p from the equation of state, 
compute g AB from (49), and t] A b from }zab and g AB ', 



and h ence o btain 7r A ' B and J 1 ' 2 . We then compute 7r a 5 



from ( 75|76 ), using v % and the equation of state. Hence 
we finally recompute 7r^v J and p — ir. We then have four 
residuals giving the discr epancy (four numbers) between 
our original guesses (192) and the recomputed values, as 
a function of the four initial guesses. We can now use any 
standard root-finding method, such as a Newton solver, 
to find the solution (and hence the correct primitive val- 
ues) to desired accuracy. (Note: this will converge only 
with a good initial guess, and a solution may not exist or 
be unique in general.) 

Note that in the fluid limit we only need to guess p, and 
our scheme then reduces to the standard conserved-to- 
primitive conversion, requiring a root find in one variable 

El- 



VI. NUMERICAL TESTS 
A. Description of the code 

The computer code used to produce the following re- 
sults uses planar symmetry; all of the variables in the 
code are 3-dimensional, but the system is only evolved 
in one or two dimensions. The numerical methods em- 
ployed are those in [22]. Briefly, the code uses a HRSC 
method with a third-order Runge-Kutta time evolution. 
In the reconstruction, standard slope limiting techniques, 
applied to the primitive variables are used - all results 
shown used van Leer's MC limiter ([23]). The HLL ap- 
proximate Riemann solver ([24]) is used to calculate the 
fluxes. The code can be run using either the relativistic 
or Newtonian set of governing equations. 

The HLL flux is 



ftqf^+ftqfHWtqf-x-qf) 



(195) 



a complete set of primitive matter variables, but these 



where q^_ x and qf are the right and left reconstructed 
vectors of conserved variables for the (i — l) th and i th 
cells, respectively, and Ahll is an estimate of the abso- 
lute value of the largest coordinate characteristic speed. 
We set this either to max |A| at one point (using the nu- 
merical calculation outlined in Sec. IV H), to max |A| over 
the whole grid, or to a constant (for example, Ahll = 1 
in highly relativistic situations and assuming the matter 
evolution is causal). 

In two dimensions standard directional splitting tech- 
niques are used. Specifically, on our (logically) Cartesian 
grid we compute the app ropri ate one dimensional fluxes 
T % required by equation (184) by sweeping through the 
grid lines one dimension at a time. The update terms 
are accumulated and applied simultaneously to minimize 
symmetry errors caused by the splitting. 

We briefly note that the performance of the code has 
been compared to a relativistic hydrodynamics code by 
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reducing the elasticity code explicitly to the hydrody- 
namic limit. As none of our codes have been optimized 
for performance, any comparisons will be approximate. 
Nevertheless, as the elasticity code is approximately 8 
times slower than the hydrodynamic code on the same 
problem (equivalently, the run-time of the hydrodynamic 
code is approximately 12% of the run-time of the elastic 
code) we see that performance will likely be an issue in 
realistic simulations. 



riD versus n^ 



As discussed in Sec. VB[ the particle number density 
n can either be obtained from the conserved variable ip A i 
(or its inverse F z a in the mixed framework), or from the 
conserved variable D. If we evolve D as a dynamical 
variable and use rip to represent the primitive variable, 
we have one more variable than if we use n^. However, 
we have shown in Sec. |IID| the evolution equations for 
no and n^jj are equivalent even if the constraints are not 
obeyed, and so we expect that both formulations have 
identical stability properties. 

In fact, when these two evolutions are compared, the 
RMS relative error in the resulting data is small; we ex- 
pect that this is finite-differencing error, as it converges 
away between first and second order. Hoewever, when 
n^p and rip are compared for a single evolution where rip 
is dynamical, the difference is of the order of round-off er- 
ror rather than finite-differencing error; we suspect that 
this is an artifact of planar symmetry. 



C. ip versus F 

A mixed framework using the inverse F A i of the con- 
figuration gradient ip l a is outlined in Appendix [d] For 
constraint satisfying initial data the results of the two 
frameworks should be the same. We have implemented 
both frameworks numerically and compared them. We 
find that the difference is on the order of the finite dif- 
ferencing error for constraint satisfying initial data, as 
expected. 

Some of the tests given in [12] and [13] do not satisfy 
the constraints - namely the second test in [12] and the 
fifth test in [13] . As the evolution of such data depends 
on the choice of constraint addition to the equations, 
we would expect it to depend on the framework used. 
Our numerical results obtained in the "mixed" frame- 
work (presented in Appendix [D|, which is that used in 
[12j[T3], matches their numerical results for all tests. Our 
results using the Eulerian framework (presented in Sec.[TT| 
and also used by [14] [21]) match only for the physical 
tests, where the initial data obey all constraints. 



D. Newtonian Riemann tests 

To validate our Newtonian code, and the Newtonian 
limit of our relativistic code, we have compared our re- 
sults to two previously published studies ([12] and [13]). 
These results use the Newtonian theory and the mixed 
framework outlined in Appendix [D] 

Broadly the results obtained from our codes matched 
those shown in [12] and [13]. As an example, we show 
the results for the first test of [12] in Figs. [6JJ9| using 
the results of the Newtonian code. The precise initial 
data used is outlined in Appendix [j] We see the seven 
waves expected for this solution; three left travelling rar- 
efactions (the second is very small), a contact, two right 
travelling rarefactions (again the second is very small), 
and a fast shock. For clarity, the wave structure of the 
exact solution is shown in Fig. [5] All waves are captured 
with only minor under /over shoots, and the numerical 
solutions converge to the exact solution [25] with resolu- 
tion, as seen by comparing Figs. [6] and [7] with Figs. [8] and 

m 

Similar results are seen for all comparison tests. How- 
ever, not all of the tests run robustly for all numerical 
methods possible within our code. An example is the 
sonic point test problem outlined in [13] (see in partic- 
ular Figs. 5-7 there). At the contact discontinuity there 
is an unphysical "dip" in the density and a correspond- 
ing "jump" in the internal energy. This is the classi- 
cal "wall-heating" effect seen by most numerical meth- 
ods when strong rarefactions separate (e.g., on reflections 
from walls or the origin in spherical symmetry - see [26] 
for the classical case and [27] for a brief discussion of 
the relativistic case). In our case these artefacts lead, 
for certain choices of numerical parameters, to numeri- 
cal results that are unphysical. This typically manifests 
itself by an imaginary characteristic speed, usually as 
the squared sound speed becomes negative. Variants of 
the code which rely on calculations of the characteristic 
information immediately fail. This problem will only af- 
fect some numerical methods at low accuracy in certain, 
somewhat artificial, situations, so is unlikely to cause 
problems in realistic situations. 

Finally, we note that a direct and comprehensive com- 
parison to the results of [13 is complicated by two is- 
sues. Firstly the units for the entropy appear inconsistent 
there, as detailed in Appendix [j] Secondly we do not find 
agreement in the comparison of the pressure tensor pij 
(denoted a there). As all other values and wave speeds 
match up well, and we have comprehensive quantitative 
agreement with the results of [12 , we believe our results 
to be correct. 



E. Newtonian limit vs Newtonian code 

The code (both relativistic and Newtonian) uses geo- 
metric units where the speed of light is one. In particular, 
all velocities are of the form v = v/c where v is a dimen- 
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FIG. 5. The density and specific internal energy in the first 
Newtonian Riemann test from [12] (from now on BDRT1), 
illustrating the seven waves possible in elastic matter. This 
is the exact solution, illustrating the wave structure in detail. 
The rarefactions - the 1, 2, 3, 5 and 6- waves - are given by the 
blue dashed lines and are shaded beneath to show the width 
of the fan. The contact - the linear 4- wave - is given by the 
dotted black line. The shock - the 7- wave - is given by the 
thick dash-dotted black line. It is clear that resolving some of 
the rarefaction waves will be difficult at moderate resolution. 



sionless velocity, v its value in conventional units and c 
the value of the speed of light in the same units. All 
parameters in the equation of state, such as e and c 2 s1 
are treated analogously. There is no need to rescale rest 
mass and length, as long as units are used consistently. 

Changing c while keeping v etc. fixed is a trivial scale 
invariance of the Newtonian equations and their solution, 
but in the relativistic equations decreasing c with v etc. 
fixed makes the same test problem more relativistic. We 
can use this to obtain an insight into the effects of (spe- 
cial) relativity, and to verify that our relativistic code has 
the correct Newtonian limit as c —> oo. 



In Fig. 10 we show the results from the relativistic code 
run with a small range of values for c. Only relatively 



FIG. 6. Numerical solution of the BDRT1 test. We show the 
results of our Newtonian code using 200 points (only 100 are 
plotted for clarity) , with the exact solution given by the solid 
line. Density, specific internal energy and normal velocity are 
shown. Only minor under /over shoots are visible. 



small values of c are shown - 3, 5, 10 and 20 km s~ , com- 
pared to a typical velocity in the (Newtonian) Riemann 
problem of 1 km s~ - as for sufficiently large values of c 
the results are visually indistinguishable. We see that the 
results from the relativistic code are qualitatively similar, 
in terms of wave structure and accuracy, and approach 
the Newtonian results in the limit c — >• oo. 



F. Relativistic Riemann tests 

In the genuinely relativistic limit we have tested our 
code against exact solutions constructed by solving a pre- 
determined wave structure. The explicit procedure is de- 
tailed in Appendix [I] and follows the method used in the 
Newtonian case outlined in [12], without constructing a 
full Riemann problem solver. 

We have verified that the code behaves correctly for 
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FIG. 7. Numerical solution of the BDRT1 test. We show the 
results of our Newtonian code using 200 points (only 100 are 
plotted for clarity) , with the exact solution given by the solid 
line. Components of the configuration gradient are shown. 
Only minor under/over shoots are visible. 



single shocks and rarefactions in the relativistic limit, 
and for some invented initial data sets that test a range 
of wave structures. As an example, we show in Figs. [T2] - 
15 the results for a four wave problem. For clarity, the 
wave structure of the exact solution is shown in Fig. [TT] 
There are two left-going rarefactions (1 and 2-waves), one 
right-going rarefaction (a 6- wave) and a right going shock 
(7- wave). The central three waves - the nonlinear 3 and 
5-waves and the contact - are all trivial. We note that 
some of the quantities change so rapidly across some rar- 
efaction waves that they are only visually distinguishable 
from shocks at high magnification. 

Even with the violent behaviour displayed across some 
waves in this four wave test, we find our code match- 
ing the exact solution well, with no unphysical oscilla- 
tions and only minor under and overshoots that converge 
away with resolution. There are the expected minor os- 
cillations near the trivial waves, most noticeable near the 



FIG. 8. Numerical solution of the BDRT1 test. We show the 
results of our Newtonian code using 1000 points (only 100 are 
plotted for clarity) , with the exact solution given by the solid 
line. Density, specific internal energy and normal velocity are 
shown. Only minor under/over shoots are visible. Comparing 
against the results in Fig. [6] we see convergence to the correct 
weak solution. 



contact, but again these converge with resolution. 



G. Two-dimensional Riemann tests 

The constraints are trivial if all variables depend only 
on one coordinate, for example when a Riemann problem 
is aligned with the numerical grid. As a first test of 
the behaviour in three dimensions and the role of the 
constraints, we have solved Riemann problems also at an 
arbitrary angle to a two-dimensional Cartesian grid. A 
method for carrying out such 2D simulations efficiently 
is described in Appendix |F| 

We put the initial discontinuity along lines x/y = 
(our ID tests), 1,1/2 and 1/5, and use a cut along the x 
axis as an approximation to a line normal to the initial 
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FIG. 9. Numerical solution of the BDRT1 test. We show 
the results of our Newtonian code using 1000 points (only 
100 are plotted for clarity), with the exact solution given by 
the solid line. Components of the configuration gradient are 
shown. Only minor under/over shoots are visible. Comparing 
against the results in Fig. [7] we see convergence to the correct 
weak solution. 



FIG. 10. BDRT1 again, but now run using the relativistic 
code with various values of c, and compared against the New- 
tonian results. The results for p, v x (appropriately scaled by 
c) are representative of the behaviour of all quantities. We see 
that as c increases the Newtonian limit is approached. 10000 
points were used in each case and only 100 plotted for clarity. 



H. Two-dimensional Rotor tests 



discontinuity. We compare this cut, suitably foreshort- 
ened, against the exact solution. 



We have not implemented the "hyperbolicity fix" con- 
straint additions for either the kinematic or dynamical 
evolution equations. In ID the equations are symmet- 
ric hyperbolic anyway, as there are no constraints then, 
but in 2D our equations are not even strongly hyperbolic. 
The error at the same time is somewhat larger in 2D than 
in ID, see Fig.[l6j but there is no sign of numerical insta- 
bility in 2D. We have no explanation for this, but expect 
that constraint addition will be necessary in other tests. 



To study a genuinely two-dimensional problem we con- 
sider a test suggested by [32] . The initial data, detailed in 
Appendix [jj represents an elastic rotor problem, where 
an inner rotating bearing is instantaneously welded to 
the non-rotating exterior, causing the rotor to slow and 
propagating elastic waves through the material. In all 
cases the rotor has coordinate radius 0.1, whilst the ex- 
terior is at rest. In all numerical experiments shown here 
400 2 points were used. 

Results for the Newtonian case are shown are represen- 
tative coordinate times are shown in figures [l7| and 18 
These should be compared to the results shown by Dumb- 
ser et al. in [32]. The results in the literature use a con- 
siderably more accurate numerical method, which is both 
higher order and uses finite elements better adapted to 
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FIG. 11. The density and shear scalar for the relativistic 4- 
wave test. This is the exact solution, illustrating the wave 
structure in detail. The rarefactions - the 1, 2, and 6-waves 
- are given by the blue dashed lines and are shaded beneath 
to show the width of the fan. The very narrow 2 and 6-waves 
are shown in detail in the insets. The contact - the linear 
4- wave - is trivial. The shock - the 7- wave - is given by the 
thick dash-dotted black line. It is clear that resolving some of 
the rarefaction waves will be difficult at moderate resolution. 



FIG. 12. Numerical solution of the relativistic 4- wave test. 
Density, specific internal energy and normal velocity are 
shown. The 4 wave structure (two left-going rarefactions, 
one right-going rarefaction and one right going shock) is diffi- 
cult to see in these variables. The solution is computed using 
200 points but only 100 are plotted for clarity. We see that 
all waves are captured well and with only minor under/over 
shoots, most visible for the second rarefaction wave in quan- 
tities such as e. 



the symmetry of the problem. Despite this, we see qual- 
itative agreement in the waves emitted during the evolu- 
tion of the problem. 

Results for the relativistic case are shown are represen- 



tative coordinate times are shown in figures 19 and 20 
Again we see qualitative agreement in the emitted wave 
structure, despite the differences in the models. 



VII. CONCLUSIONS 

We have presented a framework that can be used for 
simulating nonlinear elasticity in numerical relativity, 
and checked its viability in Riemann tests. The frame- 



work can be directly related to existing approaches and is 
a first step towards the simulation of neutron star crusts. 

Our numerical simulations show that the results from 
the Newtonian limit of the equations match those in the 
literature, and that the Newtonian limit of the relativistic 
code also match the results from the Newtonian litera- 
ture. 

The equations in first order form consist of three 
groups: evolution equations for a configuration gradient 
ip A i, auxiliary constraints ip A [ij] = for this variable 
due to the fact that ip A i = d\ A /dx l for an implicit un- 
derlying configuration x A , and energy-momentum con- 
servation laws. 

The first two groups are purely kinematical in the sense 
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FIG. 13. Numerical solution of the relativist ic 4- wave test. 
Components of the configuration gradient are shown. The 
4 wave structure (two left-going rarefactions, one right-going 
rarefaction and one right going shock) is most clearly seen in 
the plot of ip Y x . The solution is computed using 200 points 
but only 100 are plotted for clarity. We see that all waves are 
captured well and with only minor under /over shoots. 



that they are independent of the geometry of both space- 
time and matter space, and hence are the same in New- 
tonian and relativistic elasticity. However, from a space- 
time point both the evolution equations and constraints 
naturally arise as components of a single spacetime con- 
straint ip A [ a ,b] — 0- (I n f ac t, without the benefit of this 
point of view, some of the constraints seem to have been 
systematically overlooked in the Newtonian literature, 
giving rise to the numerical solution of unphysical Rie- 
mann problems in [l2j[T3].) 

Energy-momentum conservation is due to time and 
space translation invariance, and this fixes their correct 
weak form [11] . The weak form of the kinematical equa- 
tions appears to have been assumed ad hoc in the Newto- 
nian literature. Here we have rigorously derived it from 
the absence of dislocations in the elastic matter. 



FIG. 14. The relativistic 4-wave test again, but now using 
1000 points (only 100 are plotted for clarity). We see that 
all waves are captured well and with only minor under/over 
shoots, and comparing to Figure^] we see the expected con- 
vergence. 



There are two rather different frameworks in the New- 
tonian literature. One of these [14] [T5l [21] is fully Eule- 
rian and arises naturally as the Newtonian limit of our 
relativistic framework. The other [T6UT8] mixes Eulerian 
and Lagrangian points of view and gives rise to more 
complicated evolution equations. For completeness, we 
have proved that the two frameworks are equivalent in 
their weak form, and hence that the weak form of the 
second framework is also correct. This is borne out by 
our numerical tests, which agree for both frameworks (if 
the initial data obey the constraints). 

The dynamical equations of our framework can be re- 
lated to the standard Valencia formalism for relativistic 
hydrodynamics. Although, as noted in section [H] the 
fluid limit is singular, the system presented takes the form 
of the Valencia equations with additional terms. We also 
note that steps within the numerical code, such as the 
conversion to primitive variables outlined in Sec. VB[ 
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FIG. 15. The relativist ic 4- wave test again, but now using 
1000 points (only 100 are plotted for clarity). We see that 
all waves are captured well and with only minor under/over 
shoots, and comparing to Figure [l2| we see the expected con- 
vergence. 



tend towards standard algorithms in the fluid limit. 

Using the methods of [10 J, we have shown that our 
framework can be made symmetric hyperbolic, at least 
in a neighbourhood of the unsheared state of the matter, 
if certain linear combinations of the auxiliary constraints 
are added as source terms to the conservation laws. We 
have also shown that if constraints are added only to the 
kinematic evolution equations, bringing them into the 
form ([l7| (the "hyperbolicity fix"), but not to V ' a T ab = 
0, the resulting first-order system is strongly hyperbolic 
but not symmetric hyperbolic. 

The latter is precisely the situation in the Newtonian 
literature, and so the Newtonian limit of our result shows 
that the equations given there [HI [18] are only strongly 
hyperbolic but could be made symmetric hyperbolic by 
a simple constraint addition. 

There remain two outstanding issues before this frame- 
work can be used in a fully nonlinear GR simulation of 



FIG. 16. Results for the BDRT1 Riemann test, calculated on 
a two-dimensional grid for three different angles between the 
initial discontinuity and the grid, each at two resolutions. The 
initial discontinuity was placed on the line given above each 
plot. In order to compare the results to the exact Riemann 
solver presented in [12 , a slice through the two-dimensional 
grid is taken along the x axis (as an approximation to a line 
normal to the waves) , and x is scaled to correspond to distance 
perpendicular to the initial discontinuity. The spatial resolu- 
tion is independent of the angle of the initial discontinuity, 
and the snapshot is always taken at the same time, for all an- 
gles. The relativistic code is used in the Newtonian limit (as 
the exact solution is Newtonian). The high-resolution results 
were produced using Ax = Ay = 0.001. The low-resolution 
version was produced using Ax = Ay = 0.01. (For clarity, 
only 1 in 2 or 1 in 20 points are plotted for the x/y — and 
x/y — 1/5 cases, while 1 in 3 or 1 in 28 points are plotted for 
x/y = 1.) All three evolutions look similar, with the results 
approaching the exact solution as the resolution is increased; 
the only notable feature is that the left-most rarefaction wave 
does not appear to be well captured by the evolution with a 
slope of 1. 
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FIG. 17. Newtonian rotor test, following [32]. An initially 
rotating central cylinder is slowed by the interaction with the 
exterior, which is initially at rest. These figures show pv x 
at coordinate times t = 0.02,0.05,0.1 and 0.15. The results 
qualitatively match those in Figure 24 of [32] , 



FIG. 19. Relativistic rotor test, to be compared with the 
Newtonian results in figure |17[ An initially rotating central 
cylinder is slowed by the interaction with the exterior, which is 
initially at rest. These figures show pv x at coordinate times 
t — 0.02,0.05,0.1 and 0.2. The emitted waves are qualita- 
tively similar to the Newtonian results. 
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FIG. 18. Newtonian rotor test, following [32]. An initially 
rotating central cylinder is slowed by the interaction with the 
exterior, which is initially at rest. These figures show pF y y 
at coordinate times t = 0.02,0.05,0.1 and 0.15. The results 
qualitatively match those in Figure 25 of [32] . 



FIG. 20. Relativistic rotor test, to be compared with the 
Newtonian results in figure |18| An initially rotating central 
cylinder is slowed by the interaction with the exterior, which 
is initially at rest. These figures show pF y y at coordinate 
times t = 0.02, 0.05, 0.1 and 0.2. The emitted waves are qual- 
itatively similar to the Newtonian results. 
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a neutron star. The first is the issue of the integrabil- 
ity constraints in higher dimensional simulations. In the 
Newtonian literature it is clear that the hyperbolicity 
fix included here is required to obtain stable evolutions. 
However, there is no agreement as to the impact of con- 
straint violations on the accuracy of the simulations. In 
analogy with MHD simulations where the V • B = 
constraint is crucial for accuracy, we might expect that 
methods for reducing constraint violations (such as the 
parabolic damping term used by [28 - similar to the 
Powell method for MHD), or alternatively a discretisa- 
tion that maintains a discrete version of the constraints 
along the lines of Appendix [EJ will be important. 

Secondly, to be useful for simulating a neutron star, 
we must couple the elastic crust to the fluid interior. A 
framework for the nonlinear simulation of multiple mat- 
ter models separated by sharp interfaces in GR was stud- 
ied in [22] , but only for fluid- fluid interactions. This 
model built on standard Newtonian methods which have 
themselves been extended to deal with solid-fluid inter- 
actions; we expect that these methods will extend to rel- 
ativity as well. 
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Appendix A: 3+1 split of spacetime 

For reference, we assemble some standard formulas. In 
3+1 numerical relativity, the spacetime metric g a b is split 
into a spatial metric 7^ with inverse 7* J , a lapse a and 
shift ft, as 



#00 = -a 2 + ftft, goi = ft, 



7ij, 



(Al) 



where we define indices on ft to be moved implicitly with 
jij. The (absolute value of the) determinant of the 4- 
metric is given by 



9x= a lx- 



(A2) 



and hence the volume forms on M 3 and M 4 are related 

by 



Stijk „ Ajk 

e J = ae J . 



(A3) 



The inverse 4-metric is 

g 00 = -oT 2 , g 0i = oT 2 ft, g ij = j ij - oT 2 f3 i f3 J . 



The covector normal to the surfaces of constant t has 
components 

no = —a, rii = 0, (A5) 

and hence 

n° = a~\ n i = -a~ 1 l3 i . (A6) 

Hence the projector into the surfaces of constant t 

lab := 9ab + n a n b (A7) 

has components 

Too = ftft, 7oi = ft, lij = liji (A8) 

and 

7 00 = 0, 7° i= 0, 7^=7^ (A9) 

We define the convective derivative to be the derivative 
along the 4- velocity, 



n d d , d 

ox a at ox 1 



(A10) 



The factor of proportionality is given by the normalisa- 
tion condition 



u a u b g ah = -1. 



(All) 



We find 



where 



u a = {u\ u l ) = aT 1 W(l, v% (A12) 

u a = (u u Ui) = W(-a + Vjf3 J ,Vi), (A13) 



v l : = av % - ft, 

W := (1- vy) 



-1/2 



(A14) 
(A15) 



and where we define the indices on v % to be moved im- 
plicitly with jij. The scalar 



u a n n = W 



(A16) 



gives the Lorentz factor of the relative velocity between 
the matter and the time slices. 



Appendix B: Definitions of hyperbolicity 

We summarise some standard definitions [lOj [29] in 
our notation. Let w a be a vector of variables obeying 
the system of first-order partial differential equations 



P a p c wf c + l.o. = 0, 



(Bl) 



(A4) 



where l.o. stands for lower order terms. Obviously the 
index a labelling the equations needs to take as many 
values as the index j3 labelling the variables. 



24 



Assume, however, that a is an index of the same type 
as j3 and that P a (3 c = P/3a c - Then we have a conserved 
current (up to lower order terms) in the sense that 



J c 



Lo., J c := P a(3 c w a w f3 . 



(B2) 



If furthermore there exists a covector t c with the property 
that 



E(w,w) := t c J c = t c P ai3 c w a wt 



(B3) 



is positive definite, called a subcharacterisic vector, then 
the system is called symmetric hyperbolic. (In a rela- 
tivist ic context we expect t c to be timelike.) E allows 
us to estimate an L 2 norm called an energy norm of the 
solution in terms of the initial data and boundary data. 
A characteristic direction is a covector k c such that 



det k c P a(3 c = 



(B4) 



and the corresponding characteristic variable w a is the 
non-zero vector obeying 



c n „P 



k c Pa? c w p = 



(B5) 



This means that a plane wave with amplitude w^ and 
wave number k c is a solution of the principal part. For a 
causal system in relativity, influence cannot travel faster 
than light, and so k c must be spacelike or null. 
For a second-order system 



Pa0 Cd ^ cd + \.O. = 



(B6) 



the equivalent definition of a characteristic direction and 
variable is 



Cdnnf* 



k c k d P a p ca wP = 0, 



(B7) 



and it has the same interpretation as a plane wave solu- 
tion of the principal part. 

It is often useful to decompose the characteristic equa- 
tion with respect to a preferred hypersurface. Let 



k, n 



Xn„ 



(B8) 



where n a is a unit timelike covector and e a a unit space- 
like covector normal to n a . A is called the characteristic 
velocity (relative to n a ) of the characteristic variable w a . 
k a is normal to the characteristic plane spanned by the 
vectors 



Xe a 



(B9) 



where s a is any vector normal to both n a and e a (so 
that k a v a = 0). The relative speed between n a and v a 
(calculated from n a v a /\n\\v\) is ^/A 2 + s a s a > A. The 
disturbance itself moves along n a + Ae a , that is in the 
direction e a with speed A as measured by n a observers. 
One natural choice of n a is the unit normal to the surfaces 
of constant time t, and the resulting values of A are used 
in the numerical scheme. By contrast, choosing n a = u a 



gives the speed of the disturbances relative to the matter, 
which are simpler to compute. 

To make contact with non-relativistic concepts of hy- 
perbolicity, we rewrite the first order characteristic equa- 
tion (lB4|) as 



P e w = Xw, V e :=(n a P a )-\e b P b 



(BIO) 



where we have not written the Greek indices for simplic- 
ity. (If n a is subcharacteristic, n a P a is positive definite 
and so has an inverse.) The system is then called weakly 
hyperbolic with respect to the time direction n a if V e 
has real eigenvalues A for all unit vectors e a normal to 
n a . It is called strongly hyperbolic if furthermore V e 
has a basis of real eigenvectors that depends continu- 
ously on e a . It is called symmetric hyperbolic if V e is 
symmetric. As a real symmetric matrix is always diago- 
nalisable with real eigenvalues, symmetric hyperbolicity 
implies strong hyperbolicity. More generally, the system 
is called symmetric hyperbolic, or symmetrisable, if there 
exists a symmetriser, a positive definite symmetric ma- 
trix % independent of e a such that V e l~L is symmetric. In 
this case E = % a ^w a w^ . 



Appendix C: The Newtonian limit 

We obtain the limit of Newtonian motion in the ab- 
sence of gravity in two steps. In the first step, we let the 
spacetime go to Minkowski spacetime in adapted coordi- 
nates, 



ds 2 — —dt 2 + jij dx % dx\ 



(Cl) 



where 7^ is flat and independent of £, but x l could still 
be curvilinear coordinates. Hence 

v i =v\ (C2) 

and the advection equation ([29| becomes 

(d t + v i d i )k AB = 0. (C3) 

In the second step, we use dimensional analysis of the 
special relativistic equations of motion to insert a param- 
eter c representing the speed of light, as follows: 



n, 

— 2 —2 —2 

C 6, C p, C TTij, 

—3 —4 

C 7T0i, C 7Too, 



for the primitive variables, and 
D, c" 1 ^, c~ 2 t, 



(C4) 
(C5) 
(C6) 
(C7) 

(C8) 



c-^iDf, c- 2 F{Sj)\ c- 3 T(t)\ (C9) 

for the conserved variables. We then take the limit c — » 
00 of the relevant equations for Minkowski spacetime. In 



25 



this limit, 




W = l, 


(CIO) 


u a = n a , 


(Cll) 


h a b = lab, 


(C12) 


^ A t = o, 


(C13) 


TT = l lJ 7Tij = 0, 


(C14) 


D = n, 


(C15) 


Si = nvi, 


(C16) 


r = n(v 2 /2 + e), 


(C17) 


F{D) 1 = nv\ 


(C18) 


F(Sj) 1 = nvjV 1 +pS Z j + 7r*j, 


(C19) 



T(ry = n(v 2 /2 + e)v l + pu* + tt> j , (C20) 

where v l and 7r^ are now the Newtonian velocity and 
stress tensor, and their indices are moved implicitly with 
the metric jij of Euclidean space. Instead of requiring 
p, /i and fi as functions of h (the relativistic enthalphy, 
which includes the rest mass energy) and n, we need them 
as functions of e and n. The reconstruction of n, vi and 
e from D, 5^ and r becomes explicit for the equations of 
state we consider. 



Appendix D: The mixed framework 

Variables In the alternative Newtonian framework 
of [TTJ [18], the deformation is given by a map from a 
3-dimensional matter space and time to 3-dimensional 
space 



with derivatives 



F: 



dx l 






dx % 

~dt 



(Dl) 
(D2) 



(D3) 



where F l a is the 3x3 matrix inverse of i\) A %. We shall call 
this the mixed framework, as the dependent variables are 
Lagrangian, but the independent ones are still Eulerian. 
(A purely Lagrangian framework also exists, but is not 
relevant for us because we are interested in finite volume 
methods for weak solutions.) 

For the purpose of a systematic derivation of the kine- 
matic equations, and a comparison with the Eulerian 
framework, we add a time coordinate r to matter space, 
which now has coordinates ^ a — (r, ^ A ). To make this 
extension trivial, we then fix r = t. Note that 



d__ d_ 
dr ~ dt 



dx l 



(D4) 



gives us the extended derivatives 
F a n - 



Vc 



1 A 
v l F { A 

1 0< 



(D5) 
(D6) 



which are now 4x4 matrix inverses of one another, as- 
suming (10). 

Kinematic equations We derive the evolution equa- 
tions and constraints in the mixed framework by working 
in the 4-dimensional notation at first. The integrability 
condition 



«C~V := F a [a ^ = 



(D7) 



can be written as the commutator of the vector fields d a 
and d/3 pushed forward to spacetime, 



(2)/m 
^ a/3 

We define the determinant 



F a [aF /3],a = 0- 



^ 'x£, • — a I Oabcd 



5 a ^ & F a a F b pF c 1 F d s, 



(D8) 



(D9) 



where the suffixes indicate that this depends explicitly 
on the coordinates x a and £ a . With ip the inverse of F, 
we have the variation-of-determinant rule 



SR 



Xt; 



F x £ 1/J a a . 



8F a ^ 
As 5 abed and S a ^§ are constant, we therefore have 



(D10) 



F x £,b = F x £ ip a a F a a ^. 



From (Dll) and (D7), we find that 



Ca ■= (F~ i 1 F\) a = 0, 



(Dll) 



(D12) 



while combining (D8) and (D12) we obtain 



^C\p := 2 (p- i 1 F a [a F b 0] ) a = 0. (D13) 



Developing (D5) into its first row, we find 



Fx£ = TTySijk 5 F % aF 3 bF c =• F X £ 



(D14) 



Hence we obtain the 3+1 split of the 4-dimensional con- 
straints into evolution equations and constraints: 



- {3) & AB 



Ca = (f-'F 1 ^ ,=-C a = 0, (D15) 

Cr = (F~l) t + (VF-*). =:£ = 0, (D16) 
(F^F^F^)^ =: C AB = 0, (D17) 



F-i (v j F l A - VF'a) 



is then the usual convective derivative. This extension 



=: £ l A = 0. 



(D18) 
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The remaining components 

C l At = Ca, 



(3)£l 



AS 







(D19) 
(D20) 



are redundant. Note that everything is now expressed in 
terms of F % a and v l ', and we no longer need F a A- 
The corresponding jump conditions are 



Kl F¥ A(v n -s) 



F*F"a 



F-,\v n -s) 



KlF n [AFW i B] 



F-'F n A v 11 



= 0, (D21) 

= 0, (D22) 

= 0, (D23) 

= 0, (D24) 



where the n, II notation is as in Sec. II C Note that these 



jump conditions are both more nu merous and more com- 
plicated than the jump conditions ( 23|24 ) of the Eulerian 
framework. 

Equivalence with the Eulerian framework Note that 
we have used £ and C to denote the evolution equations 
and constraints in the mixed framework, and E and C 
for the Eulerian framework. 

We have the following relations between the full and 
contracted equations for F z a. 



C A = 2^ B i C l AB , 

£ = \^ B l {v'Cb - £ 1 b) , 



(D25) 
(D26) 



and the following relations between these equations and 
those for i\) A i\ 

£'a = 6 ABC 5 ijk (i> C iv l C B jk + ^ C jE B k ) , (D27) 
>• (D28) 



knC 



AB = OABCO C 



As these relations between differential equations in- 
volve multiplication by one or more factors of ip A , which 
in general is not continuous, the corresponding jump con- 
ditions may be inequiyalent. I n part icular, it is no t clea r 
if (p22j) follows from flP24]), if ( |D2l]> foll ows from (|D23|, 
if(~ 



D23|) is equivalent to (23 



(24). However, a detailed ca' 



) or if (D 24) is equivalent to 



culation shows that all these 



relations hold. 

As an example of these calculations, consider 



K{F n A 



[S njk S A BC^ B j^ C k] 

2S ni *5 ABC ty B \^ C \\j] 



(D29) 



where in the first equality we have used the cofactor rule 
and the assumption that F is the inverse of ip, and in 
the second equality we h ave used that Sabc and S ^ k are 
continuous. From (D29) we see that (23) implies (D21) 



(as claimed above), bu t the reverse is not true. In fact, 
the right-hand side of (D29) vanishes if and only if 



[r 



a 



A h 



(D30) 



for some matter space vector a A and spatial covector fc$. 



That is why the jump condition (D23) also needs to be 
imposed. 

In the papers [H CE3 [EMS H P onl y ([D22]), (|D24|) 



and (D21) are explicitly given, but flD23) appear to have 



been overlooked. In particular, the initial data for the 
second Riemann numerical test of [12] (BDRT2) and the 
initial data for the fifth R ieman n numerical test of [T3] 
(TRT5) explicitly violate flD23| ). As noted in [28 , not 
imposing the constraints (19), or equivalently (D23), in 
full corresponds to performing surgery (of the type illus- 
trated in Fig. 0) at the discontinuity. Moreover, once the 
initial data violate the constraints, the subsequent evolu- 
tion depends on how constraints have been added to the 
evolution equations. 

Equations written in terms of the densit y We have al- 
ready noted that with ([29] ) and (37), (D16) is just particle 
number conservation ( |36[ ). Note that in weak solutions, 
we must demand that ^Jk^ is everywhere continuous, a 
property that is conserved under advection. 

F X £ can also be replaced by n in the other equations 
of the mixed framework. Defining 



Pa:= 



F l 



(D31) 



we can write (D15) and (D18) as 



{W^nf A ).=0, 

(D32) 
(W^nr A ) tt + [W^n (f A v j - Pav 1 )] tj = 0. 

(D33) 

For fixed matter space index A: these happen to be iden- 
tical with the divergence constraint and the induction 
equation for the magnetic field in the formulation [30] of 
ideal magnetohydrodyamics (MHD) in general relativity. 
Taking the Newtonian limit W = 1 , jij flat and assum- 
ing Carte sian coordinates so that y^ = 1, (36), (p32| 
and (|D33|) reduce to Eqs. (3.12), (3.26) and (3^2) of|TS"/, 



where yk^ is called p re f. Further assuming y^ = 1, 
they reduce to Eqs. (6), (9) and (lb) of [12 and Eqs. (3), 
(1) and (2) of [13] . 

In terms of f l a-> the remaining equation, Eq. (D17) can 
be written as 



hw 



Jt n P[Af B] ) 



0, 



(D34) 



where ^/k^ reappears. As we have already noted, this 
constraint is not mentioned in [12j [13j [17] [TS] . It also 
does not have an equivalent in MHD. 



Appendix E: Discrete constraint preservation 

The following class of conservative numerical schemes 
preserves a discrete version of the integrability con- 
straints. With all other numerical variables defined as 
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cell averages with, by convention, integer grid index val- 
ues, define ip A i on relevant cell faces. To initialise them 
consistently, assign values to \ A a ^ cell centres at the 
initial time (in liquid as well as solid regions). Then ini- 
tialize 

1 
Ax 

1 / A 

Ay 
1 



r 
r 



x,i+l/2,j,k 



(xt+i 



3,k 



y,i,j+l/2,k 



[XiJ + l,k 



xU 



xi hk ) 



r 



,j,k+l/2 - -^ \Xi,j,k+l ~ Xi,j,k) 



(El) 
(E2) 
(E3) 



The xfjk are usec ^ om y f° r initialisation, and are not re- 
quired afterwards. We then evolve using the conservative 
equations 



dt 



r 



x,i+l/2,j,k 



1 

Ax 



tf 



A 



Tt 



j,k) 



(E4) 



and similarly for %jj A y and ijj a z , where the numerical flux 
T A ^ k is some approximation to i(j a jV j at cell centres, 
suitably limited to enforce the TVD property. Then a 



discrete version of ip A x ,y 
vant cell edges: 



■ip 1 



0, evaluated at rele- 



1 

Ay 

1 

Ax 



(^ x,i+l/2,j+l,k ~ *P x,i+l/2,j,k) 



W 



y,i+l,j+l/2,k 



r 



y,i,j+l/2,k 



)=0, (E5) 



and similarly for the other two commutators, is obeyed at 
all times if it is obeyed initially. The fundamental idea is 
that the numerical fluxes T A are the time derivatives of 
the underlying x A •> an d hence are the same for the tp A x , 
i/j A y and ij) A z . The discrete constraints act as discrete 
integrability conditions that allow us to reconstruct the 
xfj k by summation if desired. 




FIG. 21. Example of a two-dimensional grid with shifted 
periodic boundary conditions, with n x = 8, n y = 2 and S x = 
3. The physical grid is surrounded by one ghost cell on each 
side. In reality, n x would be much larger, while n y ranges 
from 1 to a few, and S x from to a few, with no common 
factor. The placement of the left and right state of a Riemann 
problem is shown by the letters L and R and shading. Cells 
are initialised with the left or right state depending on the 
position of the cell center. The numbers 1, 2, 3, 4 identify 
four physical cells and the ghost cells they donate values to. 
The initial discontinuity is at an angle a from the y axis, with 
tana = 5 x /n y (= 3/2 in this example), and goes through the 
point x = y = 0. 



ID tests), and use the x axis as an approximation to a 
line normal to the initial discontinuity when taking a cut 
through the solution. 



Appendix G: Shear scalars 

The three eigenvalues of t] A b can be parameterised as 
{a, 6, l/(ab)}. We then find that in the unsheared state 
a = b = 1, 



Appendix F: Riemann tests on a 2-dimensional grid 

As a first test of the role of the constraints in hyper- 
bolicity, we numerically solve Riemann problems on a 
2-dimensional grid, with the initial discontinuity at an 
angle to the grid. Assume the grid consists of n x x n y 
cells, surrounded by the necessary number of ghost cells. 
After each time update, the ghost points are filled using 
periodic boundary conditions, identifying cell (i,j) with 
(i + n x ,j) in the x direction, but (i, j) with (i + S x ,j+n y ) 
in the y direction, where S x is an offset. Consistently 
with these boundary conditions, the initial discontinuity 
is then placed on a line of x/y = S x /n y (assuming that 
the grid spacing is the same in the x and y directions). 



This is illustrated in Fig. [21 

As the x and y directions are interchangeable, the 
slope 5 x /n y and its inverse pose the same Riemann test. 
(Less obviously, in our implementation those two tests 
also have roughly equal computational cost.) We choose 
n y > S x (and typically S x = 1) so that the initial discon- 
tinuity is always closer to the y axis (where it is in the 



I 1 = I 2 = 3, 


(Gl) 


I]a = I\ = II = 1% = 0, 


(G2) 


7-1 _ 7-1 _ Q rl—l 

l ,aa — 1 ,bb ~ z > i ,o6 _ ^ 


(G3) 


I,aa = ^,66 = $, ^ \ab = 4. 


(G4) 



Hence 4 (J 1 — 3) and I 2 — 3 are the same function of the 
shear up to quadratic order. This is not a bad choice of I a 
but a property of any shear invariant. It is related to the 
fact that the characteristic speeds in the unsheared state 
depend on f\ and fi only through the one comb ination 
/i + 4/2 that appears in the shear modulus (166). 



Therefore, to model linear elasticity correctly, it is suf- 
ficiently general to make the ansatz 



e(s,n,n=e(n,s) + ^^-S(n, 



n 



where the shear scalar S obeys 



s = o, 

dS 



dl 1 dP 



1 



(G5) 

(G6) 

(G7) 
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in the unsheared state I 1 = I 2 = 3, but is otherwise 
arbitrary. For any such c hoice of 5, /2(n, s) evaluates 
to the usual shear modulus (166) in the Newtonian limit, 



and the equations of motion are the same when linearised 
about the unsheared state. 

Clearly there are many possibilities of defining a shear 
scalar that obeys these conditions, but we are not aware 
of any physical reason given in the literature for why a 
specific choice should be preferred, or of values given for 
/i and J2 independently. 

An equation of state for copper in [13] uses the shear 
scalar 



1\2 



Si 



Cran - _ 



3 J 2 - (I 1 ) 
12 



(G8) 



which is homogenously quadratic in the eigenvalues of 
t) A b- In [9] the shear scalar 



<?ks := 



(n 



1\3 



lr2 



I l I 



18 



24 



(G9) 



which is cubic, is suggested for what seem to be aesthetic 
reasons. Yet another shear scalar is 



>VM 



s ab s ab 



2I 1 



where 



Sab := -Ah ah - T]ab) 



(G10) 



(Gil) 



is the "constant volume shear tensor" defined in |5]. In 
the Newtonian limit, near the unsheared state, <Svm is 
related to the Von Mises stress scalar (assuming stress 
and strain are related linearly). It gives the same values 
of fi and f 2 as <S C ran. 



Appendix H: Equations of state 

We now consider examples of equations of state of the 
form ( G5 ) . The following general expressions will be use- 
ful: 



h 



fa 



l + e + 


n n 


9 de 
on 


( dfi 


/x(n,s) 


dS 



fi)S, 



dl° 



(HI) 
(H2) 
(H3) 



In principle we can eliminate s from these two equations 
to obtain p, f\ and fi, as functions of (n, h, I 1 , J 2 ), which 
we need in the recovery of the primitive from the con- 
served variables. 

A toy relativistic EOS As a toy model for a relativis- 
tic equation of state, we take e from the commonly used 



"Gamma-law" hot equation of state, and make the shear 
modulus jl a power of the density only, namely 



e(n,s) 



m „r-i 

r-i n ' 

Kn x , 



(H4) 

(H5) 



where T, k and A are constants. This is motivated by the 
fact that in neutron star crusts \i ex n 4//3 , with the factor 
of proportionality only weakly temperature-dependent. 
The bulk modulus in neutron stars is given by the nu- 
clear interactions, while the shear modulus is provided by 
Coulomb interactions, which makes it independent and 
mu ch sm aller. Following [9], we choose S as <Sks given 
by dG9|. 



The expressions we need for the conserved to primitive 
variables conversion are then 



p(h, n, I a ) = ^^nih - 1) + ^WS, (H6) 
p(e, n, I a ) = (r - l)ne + (A - I>n A <S, (H7) 



%,n,/ a )-l + 
fi = ™ x 
h = -kt< 



r p , r-A ,_! 



r- in r-i 



24 



A-l 



24' 



nn x - l S, (H8) 
(H9) 
(H10) 



The characteristic speeds in the unsheared state are 



\rr< 



A| = 



,A-1 



1 + Te' 

r(r-i)e- 



,A-1 



1 + Te 



(Hll) 

(H12) 



Cranfield EOS The equation of state for copper used 
in [13] for Newtonian shock tube problems, translated 
into our notation, is 



e(s, n, I a ) = A(n) + B(n)K(s) + C(n)S, (H13) 

i 2 
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K 

2a 2 


A«o/ 


- 1 
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-<^' 
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= e c v - 


-i, 





C:=B t 



°\no) 



/3+4/3 



(H14) 

(H15) 
(H16) 

(H17) 
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where S is Scran given by (G8). We need the following 



forms of the equation of state: 

p(s, n, I a ) = n [nA' + jBK +{{3 + 4/3)C<S] , (H18) 



p(h,n,I a ) 



n 



j(h - 1) - jA + nA' 



7 + 1 

+(/3 + 4/3- 7 )CS 

p(e, n, 7 a ) = n 76 — 7 A + nA' 

+(/3 + 4/3- 7 )C<S 



ft(p, n, 7 a ) = 1 + ^±1 - + A - -nA' 

7 n 7 

--(/3 + 4/3- 7 )C5, 

7 

_ ^ 

C 



(H19) 
(H20) 

(H21) 

(H22) 
(H23) 



Appendix I: Constructing exact solutions 

The exact solution of the Riemann problem is a stan- 
dard test for HRSC methods. For Newtonian elasticity 
exact solvers have been constructed both by Miller [28] 
and by Barton et al. [T2] . In the relativistic case here we 
have not constructed a generic solver to compute the full 
Riemann problem solution. As noted by [12 , this can 
be extremely sensitive to initial guesses used. Instead we 
construct exact solutions by specifying the wave struc- 
ture explicitly in advance and solving across each wave. 

As summarized in [12 , with piecewise constant initial 
data the generic solution will contain seven self-similar 
waves. The central wave will be a contact discontinuity, 
and the other waves will be genuinely nonlinear. We as- 
sume that the solutions are simple shocks or rarefactions. 
We then solve across each wave in the following manner. 

Shock wave We assume that the primitive variables 
to the left of the wave, w^, are given. We then impose 
the value either of the shock speed s^ or of one compo- 
nent of the variables to the right of the wave, w#. The 
Rankine-Hugoniot conditions 



f (Wfl) - f (w L ) 



q Cp) 



[q(wfl) - q(w L )] (II) 



then form a system of nonlinear equations for the remain- 
ing components of w^ and, where necessary, for the shock 
speed s( p \ Here ^ denotes the wave number counting 
from the left. 

This problem is solved explicitly using the Matlab 
solver f solve. It is usually necessary to experiment with 
the imposed value and initial guesses in order to construct 
a solution satisfying the Lax entropy condition 



\ {p) (w L )>s^ >\( p \w R ). 



(12) 



The construction of the eigenvalues A( p ) is discussed be- 
low. 



Contact discontinuity A contact must satisfy the 
Rankine-Hugoniot conditions ( [IT] ) combined with the re- 
striction that the wave speed s matches the normal ve- 
locity on either side of the wave. Hence we can use the 
same techniques as for the shock with the value of the 
velocity imposed. 

Rarefaction wave As noted by [12] the solution across 
a rarefaction wave is given by 



dw 



M 



(w) 



r(p) (w) • V w A(p) (w) 



(13) 



Here £ = x/t is the self-similarity variable. We have that 
A^(w^) < £ where ^ labels the wave number and w^ 
is given, as above. We impose that £ < <£r = A^(w^) 
to stop the integration. In addition r^ are the right 
eigenvectors associated with the p th eigenvalue \( p \ and 
V w denotes the gradient operator with respect to the 
vector of primitive variables. 

All characteristic information (X^ p \ r^) is constructed 
from the Jacobian matrix 



J 



<9f(w) 
<9q(w) 



(V w q) X V w f. 



(14) 



As in the Newtonian case discussed in [12] we need to 
explicitly modify the calculated Jacobian to build in the 
hyperbolicity corrections as in equation (17). 



Given an explicit left state w^ the numerical solution 
is found by solving the ODE ( |l3] ) for w with initial data 
wl in A^(w^) < £ < £r. Explicitly we use the ode45 
routine with Matlab. The Jacobian J is constructed us- 
ing explicit finite differencing by varying each component 
of w by a small value h. Standard Matlab routines were 
used to construct and sort the characteristic information. 
The gradient V W A^ P ^ (w) was also constructed using ex- 
plicit finite differencing. In all cases 6 th order finite dif- 
ferencing combined with Richardson extrapolation was 
used to ensure sufficient accuracy. 

There are two potential problems with this construc- 
tion. First, as noted by [12], we have no guarantee that 
equation (|l3| has a unique solution. This would imply 
that the true solution is a compound wave, and breaks 
the assumptions made here. Second, the numerical con- 
struction of the characteristic information is extremely 
sensitive when the eigenvalues are close to each other. 
This appears to be the case for the problems and equa- 
tions of state considered below, and means that for the 
slower 3 and 5 waves next to the contact we are forced 
to construct very small rarefaction fans. 

In principle there is no reason why the procedure above 
could not be extended to construct a full Riemann solver. 
However, such a solver would have little practical utility, 
even if it could be made generic and robust. Numerical 
experiments have shown that it is faster to compute an 
approximate solution using 800 grid cells than it is to con- 
struct one exact solution with a pre-specified wave struc- 
ture. Even allowing for the massive speed improvements 
possible within our current exact solver, it is clearly im- 
practical for use within an evolution code. 
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Appendix J: Initial data for numerical tests 

We used several sets of initial data that were defined 
in published papers; this was done to ensure that our 
code agreed with Newtonian results produced previously 
[T2] [13] . Because both papers chose entropy, 5, as a 
primitive variable, instead of the pressure, p, we list the 
initial entropy value here, and calculate the pressure from 
the entropy when the system is initialized. 

For the following sets of initial data, the spacetime 
metric is the Minkowski metric, and the matter- space 
metric is the Euclidean metric in Euclidean coordinates 
normalized with the initial density of the elastic medium, 
no; we note that while we must convert units of ve- 
locity to geometrized units, we do not need to convert 
units of density or of length, as long as we are con- 
sistent throughout the code. For this paper the value 
n = 8.93 g/cm 3 was used for the BDRT tests (from [12J) 
and no = 8.9 g/cm 3 was used for the TRT tests (from 
[13]). In addition to this, for each of these situations, the 
Cranfield equation of state, described in Appendix [H[ 
was used. For comparison purposes, the velocities in this 
section are taken to be in km s _1 , while the entropy is in 
kJg^K" 1 . 

BDRT1 This is the same as Testcase 1 in [12 . It 
allows us to examine the entire seven-wave structure of 
the solution. Using the Cranfield EOS above, the solu- 
tion consists of three left-travelling rarefaction waves, a 
right-travelling contact, two right- travelling rarefactions, 
and a right-travelling shock wave. The initial data is pre- 
sented for the state vector w = (y l , F l a-, s) in the mixed 
framework given in Appendix [DJ and all other quantities 
are derived from them: 



w L 



w# 




,0.001 



(Jl) 



Results are shown at coordinate time t = 0.06. 

4-wave relativistic solution We constructed a range 
of relativistic solutions, mostly consisting of a single 
shock or rarefaction, using the technique outlined in Ap- 



pendix [T] The toy relativistic equation of state given 
above is used, with parameters T = 5/3, A = 4/3, and 
k = 1/2. In particular, we present a solution with four 
nonlinear waves. The two left-going waves (1- and 2- 
waves) are rarefactions. The contact is trivial, as are 
the central (3- and 5-waves) nonlinear waves. The slower 
right-going wave (a 6-wave) is a rarefaction, and the fast 
right-going 7-wave is a shock. The initial data is pre- 
sented for the state vector w = (v\ifj A i,p), truncated to 
6 significant figures, and all other quantities are derived 
from them: 



w L 



Wi? 





,1.86054 



(J2) 



0.469381 \ 
0.0332532 , 
0.349709 J 



0.764910 0\ 
-0.541672 1 ,0.450123 
0.369075 1/ 

(J3) 



Results are shown at coordinate time t = 0.25. 

In addition to Riemann problem style tests we consider 
a genuinely two-dimensional rotor test. The Newtonian 
rotor test was suggested by [32] , where the evolution was 
shown using a high-order finite element technique. The 
domain is cylindrical, of total radius 0.5. The material 
is initially at rest except in the rotor, represented by a 
cylinder of radius 0.1, within which it rotates with angu- 
lar velocity uj = 10. The material is not deformed (i.e., 
F l a is the unit matrix) nor hot (i.e., 8 = 0). All other 
matter properties follow the Riemann tests above. That 
is, the initial density is given by no = 8.93 g/cm 3 and the 
Cranfield equation of state, described in Appendix[HJ was 
used. Here, as we have used a Cartesian grid, we have 
simulated the full domain x,i/G [—0.5,0.5]. 

We suggest a relativistic rotor test as a direct compar- 
ison with the Newtonian version. The domain remains 
the same as the Newtonian case. The angular velocity 
is reduced to uj = 0.5. The material is initially set so 
that ip A i is the unit matrix and p = 1. As the shear 
also depends on the velocity through ^ A t, the material 
is sheared within the rotor initially, in contrast to the 
Newtonian case, but this is small. As in the Riemann 
tests above we use the toy relativistic equation of state 
given in Appendix [HJ with parameters T = 5/3, A = 4/3, 
and k = 1/2. 
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